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FOREWORD 


The computer programs described in this 
report were developed by The Analytic Sciences 
Corporation (TASC) during the period from March 3, 

1975 to November 25, 1975 under Contract No. NASl- 
13807 for the National Aeronautics and Space Admin- 
istration, Langley Research Center, Hampton, Virginia. 
The work was sponsored by the Navigation and Guidance 
Research Branch of the Flight Instrumentation Divi- 
sion as a contribution to the VTOL Automatic Landing 
Technology (VALT) Program. Dr. David R. Downing 
served as Technical Monitor for this contract. 

The author acknowledges the assistance and 
technical guidance provided by Mr. John R. Broussard, 
Mr. Paul W. Berry, and Mr. Joel M. Winett during the 
course of the computer program development . 
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ABSTRACT 


Computer programs for the design of analog 
and digital flight control systems have been devel- 
oped and are documented in this report. The pro- 
gram DIGADAPT uses Linear-Quadratic-Gaussian (LQG) 
synthesis algorithms in the design of command- 
response controllers and state estimators, and it 
applies covariance propagation analysis to the se- 
lection of sampling intervals for digital systems. 
Program SCHED executes correlation and regression 
analyses for the development of gain and trim sched- 
ules to be used in open-loop explicit-adaptive con- 
trol laws. A linear-time-varying simulation of air- 
craft motions is provided by the program TVHIS, which 
includes guidance and control logic, as well as models 
for control actuator dynamics. The programs are coded 
in Fortran and have been compiled and executed on both 
IBM and CDC computers. 
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1 . 


PROGRAM OVERVIEW 


This report describes computer programs which can 
be used in the analysis, design, and evaluation of flight 
control laws. Design techniques for analog (continuous-time) 
and digital (discrete-time) control systems are implemented 
in Program DIGADAPT using Linear-Quadrat ic-Gausi^ian (LQG) 
Synthesis. Program SCHED computes correlations between con- 
trol constants and flight conditions, indicating possible 
relationships for scheduling gain and trim settings in the 
actual system. The program also generates the gain schedul- 
ing functions by regression analysis. Program TVHIS is a 
linear-time-varying simulation of flight motions in the neigh- 
borhood of a reference flight path. It is a useful tool for 
evaluating the performance of control systems alone or in com- 
bination with predetermined guidance laws. 

Design procedures using these programs can be flexi- 
ble, as several options for estimation and control logic are 
provided. Control laws include dynamic, proportional-ir- * ' gral 
(PI), and proportional-double integral (PII) command-response 
controllers, as derived in Ref. 1 and summarized in succeeding 
chapters. Estimation logic is based upon the Kalman filter, 
and the sampling interval for digital systems is determined 
using a state covariance propagation technique. Once system 
gains have been computed at several flight conditions, the 
gain scheduling feature allows an open-loop explicit-adaptive 
control law to be defined. 

These design programs have been developed for a spe- 
cific application, the control of a tandem- rotor helicopter, 
and they can be extended to a variety of aircraft configurations, 
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including other VTOL designs and conventional transport, 
fighter, and general aviation aircraft. Application to re- 
entering spacecraft, missiles, and submarines also is straight- 
forward. In each case, details of control effectors and elas- 
tic modes must be specified by the designer. The basic routines 
and subroutines have been programmed to handle systems of arbi- 
trary dimension, go they are generally applicable to a wide va- 
riety of vehicle cQntrol..la^^ 

1.1 GENERAL SYSTEM DESCRIPTION 

DIGADAPT and SCHED are used to develop a flight con- 
trol structure and a gain scheduling table. The purpose of 
TVHIS is to evaluate and determine the success of the design 
procedure. Figure 1.1-1 gives an overview of the control flow 
in the design procedure, and Fig. 1.1-2 gives the data flow 
for this system. 

1.2 PROGRAM IMPLEMENTATION 

The design programs were developed on an IBM 370/145 
using the Fortran H compiler. They also have been converted 
to run on a CDC 6600 under the SCOPE 3.3 operating system 
using the FTN Fortran compiler (Ref. 2). A number of design 
subroutines were adapted from Ref. 3. 

A permanent disk file is used to store aerodynamic 
coefficients. This file must be established using the proce- 
dure described in Appendix A before DIGADAPT or TVHIS is execu- 
ted. 

On the CDC 6600 computer DIGADAPT requires 44,352 
decimal (126,500g) locations for execution; approximately 
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40,000g locations are required for array storage to execute a 
design procedure for a 13-state, 4-cortrol system; and approxi- 
mately 30,000g locations are used by the SCOPE operating sys- 
tem. An overlay or segment procedure has not been defined 
for the following reasons : 


• Two distinct algorithms are used to solve 
matrix Riccati equations. This assures 
convergence of solutions in minimum time 
but increases the storage requirement . 

• System overhead which W5uld be incurred with 
overlay and segment structure is eliminated; 
hence, run time for each design is reduced. 

SCHED requires 33,664 decimal (101,600g) locations for 
execution and TVHIS requires 33,024 decimal (100,400g) 
locations for execution. 


1 . 3 DOCUMENTATION 

This document is divided into seven chapters. Chapter 1 
gives an overview of the three programs, DIGADAPT, SCHED and 
TVHIS, and a discussion of program implementation and documen- 
tation. Chapter 2 describes all the inputs and outputs of 
DIGADAPT, and Chapter 3 describes in detail the significant 
subroutines of DIGADAPT. Similarly, Chapters 4 and 6 give 
the inputs and outputs of SCHED and TVHIS, respectively, while 
Chapters 5 and 7 give a detailed description of SCHED and TVHIS, 
respectively . 

Appendix A is concerned with the creation of the 
permanent disk file containing aerodynamic coefficients. It 
describes how to pick appropriate cases once the file has 
been created. Appendix B illustrates the flowcharting sym- 
bols used, and Appendix C outlines the design procedure used 
in DIGADAPT and SCHED. Appendix D describes the list of de- 
liverables. 
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2. PI GAD APT USAGE 

DIGADAPT has been written using top-down structure 
and depends on user input to determine which of the program 
chapters are to be used for a given design (see Table 2-1). 
There are seven such program chapters:* The main executive, 
system matrix computation, sampling interval determination, 
linear-optimal gain design, estimator /observer design, calcu- 
lation of eigenvalues, and time history simulations. Two 
flags control the main logic, ICHPTR(7,6) and ITYP(2). 
ICHPTR(I,J) indicates which chapters are to be executed; 

I represents the i^^ chapter and J represents the section 
within the i‘' chapter. The default values for ICHPTR(I,J) 
are all 0 (representing "off”), and it is necessary for the 
user to override these values via NAMELIST by setting 
ICHPTR(I,J) to 1 (representing "on"). ITYP is a working 
flag whose definition depends upon the procedure being exe- 
cuted. Table 2-2 is a summary of flag settings for specific 
features of DIGADAPT, 


2.1 DIGADAPT INPUT 

Input to DIGADAPT takes two forms: card input on 

Fortran logical unit 5 (TAPES) and permanent file disk inputs 
on Fortran logical unit 1 (TAPEl). The card input consists 
of a header card for the first case followed by NAMELIST in- 
put cards used to override program variable defaults set in 

♦Program chapters refer to the logical structure of the pro- 
gram DIGADAPT and should not be confused with chapters of 
this document . 
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TABLE 2-1 

PROGRAM CHAPTERS OF DIGADAPT 


1. Main Executive - Run/Case Setup 
Initialize constants 

Read in dimensions, flags and case number 

2. Full Model - Stability Derivatives, F and G (body axes) 

2.1 F and G data input 
Primary routine - AERO 
Secondary routine - FGCO^^P 

3. Sampling Interval Determination 

3.1 Sampling time calculation 
Primary routine - SID 
Secondary routine - SAMPLE 

4. Digital - Adaptive Control Law 

4.1 Q and R input 

4.2 Addition of integrator states for Proportional-plus- 
Double Integral (PII) controller 

4.3 Linear-optimal regulator gains; continuousrtime case 

4.4 Discretization of continuous-time system 

4.5 Linear-optimal regulator gains: discrete-time case 

4.6 Proportional-integral (PI) controller 
Primary subroutine - GAIN 

Secondary routines - modified ESL package 

5. Estimator Design 

5.1 Kalman filter gains - discrete-time case 
Primary routine - EOD 

Secondary routines - ESTOBS, KFD 

6. Eigenvalues - Stability Criteria For LTI Models 

6.1 Eigenvalues of the closed-loop-continuous system ACL 

6.2 Eigenvalue.s of the closed-loop-discrete system ACLD 

6.3 Eigenval.ues of the closed-loop-discrete Kalman filter 
FCL 

Primary subroutine - FREE 
Secondary routines - ATEIG, HSBG 

7. Time History Computation 

7.1 LTI time history of PII controller or dynamic 
controller 

7.2 A discrete time history of PII system will be 
performed 

7.3 A discrete time history of the PII system is computed 

7.4 An approximate discrete-time history of the PII 
system is computed 

7.5 A discrete time history of the pi and/or control-rate 
system is computed 

Primary subroutine - TIME 

Secondary subroutines - TIMHIS, DTMHIS, PIDHIS 
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TABLE 2-2 

SUMMARY OF DESIGNS 



CONTINUOUS 

TIME 

r 




DISCRETE 

TIME 





□ 

DYNAMIC ! 

CONTROLLER 

PI I ATTITUDE 
CONTROL 

PI I VELOCITY 
CONTROL 

PI ATTITUDE 
CONTROL 

PI VELOCITY 
CONTROL 

PI I ATTITUDE 
CONTROL 

PII VELOCITY 
CONTROL 

DYNAMIC CONTROLLER 
VELOCITY CONTROL 

DYNAMIC CONTROLLER 
ATTITUDE CONTROL 

TRIM VALUES 

TRIM MATRICES 
F AND G 

SAMPLING INTERVAL 

KALMAN FILTER 

SAMPLE 3 

SAMPLE 2 

SAMPLE 1 

ICHPTR(2,1) 

1 

B 

B 

B 

B 

B 

B 

B 

B 

1 

B 

H 

D 

B 

B 

B 

ICHPTR(3,1) 









mm 





1 

Ki 

Bi 

ICHPTR(4,1) 

mm 

1 

1 

i 

wm 

mm 

1 

n 

B 

■ 





D 

a 

ICHPTR(4,2) 

B 

n 

n 



m 










iJ 

ICHPTR(4,3) 

wm 

1 

1 


mm 

B 









D 

a 

ICHPTB(4,4) 




kb 

B 

wm 

1 

1 

1 






n 

a 

ICHPTR(4,5) 


■■ 


1 

1 

1 

KB 

KB 

1 






a 

D 

ICHPTR(4,6) 




1 

1 



1 

KB 






1 


ICHPTB(5,1) 













1 



1 

ICHPTR(6,1) 

ma 

mm 

BB 













1 

ICHPTR(6,2) 

mm 


— 

1* 


1* 

1* 

■a 

■B 



■ 





ICHPTR(6,3) 

B 




BK 







■ 

1 

■ 

■ 

m 

1CHPTR(7,1) 

1 

KB 

KB 


mm 







■ 


■ 

a 

1 

ICHPTR(7,2) 




J 


1 

1 


, 







1 

ICHPTR(7,3) 




wm 


1* 

1* 







■ 

IK 


ICHPTR(7,4) 




B' 


1 

1 


, 





B 

Bi 

1 

ICHPTR(7,5) 




1 

1 



1 

1 



■ 



H 

m 

ITYP(l) 

1 

3 

2 

1 

1 

3 

2 

n 




■ 



B 

a 

ITYP(2) 



2 


ni 

5 


kb 

Bi 

6 

3 

Hi 

KS 


D 

a 

NPUNCH 

B 


BKI 

KB 

Bi 

1* 

KB 

!♦ 

KB 

m 

wm 

■i 

10 





♦Optional 


BLOCK DATA. For multiple case runs , header cards for subse- 
quent cases follow the NAMELIST input . Note that the NAMELIST 
input is read only on the first case run and, therefore, if 
a multiple case run is made, the NAMELIST input values will 
be the same for all cases. In addition, if a multiple case 
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run is made, case numbers must be input in increasing order 
(i.e.. Case 47 should be input before Case 55). The disk 
file input contains 297 different cases of aerodynamic 
coefficients for the system matrix F. (See Appendix A). 

All input to DIGADAPT is read by the main executive. 
Figure 2.1-1 shows card input data for a multiple case run 
and Tables 2.1-1 to 2.1-10 give a description of all DIGADAPT 
input. Table 2.1-1 describes the header card and Tables 
2.1-2 to 2.1-10 to describe namelists NAMl to NAM9 respec- 
tively. 



Figure 2.1-1 Input on Logical Unit 5 

(TAPE5) for DIGADAPT 


reproducibility of the 


TABLE 2.1-1 


HEADER CARD INPUT TO DIGADAPT 


VARIABLE 


DESCRIPTION 


L 

M 

N 

NCASE 


NSIGMA 


NDOT 


NPUNCH 


NCURV 


Number of observations 

Number of controls 

Number of states 

Case number from 1 to 297 
to determine which set of aero- 
dynamic coefficients to use, see 
Appendix A 

Flag 

=0, if SIGMA(I), 1=1 to 16 is to 
be computed in QRCOMP 
=1, if SIGMA(I), 1=1 to 16 has 
been input via NAMELIST 

Flag 

=0, if DBDOT, DCDOT, DSDOT 
DRDOT are to be used to compute 
the diagonal values of R matrix 
=1, if diagonal values of R ma- 
trix are to be read in via 
NAMELIST 

Flag 

=0, no punched output 

-1 , punched output for SCHED 

Flag - set to zero - not used 


♦I = represents a right Justified integer value 


UNITS 


TYPE* 

I col. 1-5 
I coi. 6-10 
I col. 11-15 
I col. 16-20 


I col. 21-25 


I col . 26-30 


I col. 31-35 
i col. 36-40 
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TABLE 2.1-2 


NAMl INPUT 


VARIABLE 

DESCRIPTION 

UNITS 

TYPE 

DEFAULT 

lOiPTR 

7x6 array of flags used to 
turn on (~1) and off (•■0) the 
given sections of each chapter 
where I ■ chapter and J • section 

- 

I 

43 • 0 


TABLE 

2.1-3 




NAM2 

INPUT 



VARIABLE 

DESCRIPTION 

UNITS 

TYPE* 

DEFAULT 

ROLRAT 

Uaxlmura roll rate 

deg/aee 

R 

-25 

PITRAT 

Maximum pitch rate 

dcg/sec 

R 

-25 

YAWRAT 

Maximum yaw rate 

deg/aec 

R 

-25 

DIFCOL 

Maximum differential collective “ 
pitch control 

Id 

R 

-6.5 

OOLLEC 

Maximum collective • vertical 
velocity control 

la 

R 

-4.56 

CYCLIC 

Maximum cyclic - roll control 

In 

R 

-3.6 

DIFCYC 

Maximum differential cyclic “ 
yaw control 

In 

R 

-4.18 

PITANG 

Maximum pitch angle 

deg 

R 

-10 

ROLANC 

Maximum roll angle 

deg 

R 

-10 

BEDANG 

Maximum yaw angle 

deg 

R 

-10 

XI 

If ITYP(l) " 2, then weighting on 
integrators for velocity control 

deg-a«c 

R 

1 


system - X position 
If ITYP(l) “ 3, then weighting on 
integrators for attitude system ^ 
pitch integral 




T2 

If ITVPCD - 2, then weighting 
on integrators for velocity 

ft 

R 

1 


control system - y-posltlon 
If ITYP(l) - 3, then weighting on 
Integrators for attitude system — 
roll integral 

deg-aac 

R 


T3 

If ITYP(l) - 2, then weighting 
on integrators for velocity 

ft 

R 

3 


control system - z-posltion 
If ITTPCD « 3, then weighting 
on Integrators tor attitude 
system - yaw Integral 

deg-sec 

R 


T4 

If ITYP(l) ” 2, then weighting 

deg-sec 

R 

1 


on Integrators for velocity 
control system - yaw Integral 
If ITYP(l) - 3, then weighting 
on Integrators for attitude 
system - z-positloo 

ft 

R 



•R represents a real variable 
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TABLE 2.1-4 
NAM3 INPUT 


VARIABLE 

DESCRIPTION 

UNITS 

SIGMA 

16 vector - values used in com- 



putlng diagonal values of Q 
matrix (see subroutine QRCOMP) 


KFLAG 

9 vector flag to determine 
which states are assumed to be 

- 


measured In discrete form 



N> 

I 


TABLE 2.1-5 
NAM4 INPUT 


VARIABLE 

DESCRIPTION 

UNITS 

IXX 

Second moment of inertia about 
x-axis = 

slug-f t^ 

lYY 

Second moment of inertia about 
y- axis - 1^^ 

slug-ft^ 

slug-ft^ 

IZZ 

Second moment of inertia about 

z-axls = I 

zz 

IXZ 

Product of inertia = I 

xz 

slug-ft'* 


» I 


TYPE 

DEFAULT 

R 

- 

I 

9*1 
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TABLE 2 . 1-6 
NAM5 INPUT 


VARIABLE 

DESCRIPTION 

UNITS 

TYPE 

DEFAULT 

DBDOT 

Pitch control rate = 

in/sec 

R 

>2.0 

DCDOT 

Vertical velocity control 
rate - 5 ^ 

in/sec 

R 

-2.0 

DSDOT 

Roll control rate =• 6g 

in/ sec 

R 

-2.0 

DRDOT 

Yaw control rate =• 6^ 

in/sec 

R 

-2.0 

UDOT 

U acceleration * u 

ft/sec^ 

R 

-0.2 

VDOT 

V acceleration » v 

ft/sec^ 

R 

-0.2 

WDOT 

W acceleration ■ w 

ft/sec2 

R 

-0.2 

UMAX 

U velocity ■ u 

ft/sec 

R 

-1.5 

VMAX 

V velocity ■ v 

ft/sec 

R 

-0.5 

WHAX 

W velocity ■ w 

ft/sec 

R 

-0.2 


TABLE 2.1-7 
NAM6 INPUT 


VARIABLE 

DESCRIPTION '• 

UNITS 

TYPE 

DEFAULT 

Q3AR 

Noise covariance matrix dimen- 
sion (N+3) X (N+3). Note 
singly subscripted* 

- 

R 

256 * 0.0 

SVB 

13 state vector bound 

- 

R 

2. .2. ,2. .5. .5. ,5.. 
.5..5..5.,1..1..1..1. 

SIG\fAD 

Covariance matrix for process 
noise (9x9) 

- 

R 

81 ♦ 0.0 

THETAD 

Covariance matrix for observation 
noise (Lx9) - Note singly sub- 
scripted* 


R 

290 * 0.0 


♦To input a value in the (I,J)th position of the matrix, use the following algebraic conversion to 
compute appropriate subscript for singly subscripted variables. 


K = (J-1)* N + I 

where N is the number of rows of the matrix 
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TABLE 2.1-8 
NAM7 INPUT 


VARIABLE 

DESCRIPTION 

UNITS 

TYPE 

DEFAULT 

P 

(NxN) system matrix - must use 
single subscript* 

- 

R 

- 

G 

(NxMJ control matrix - must use 
single subscript* 

- 

R 

• 

Q 

(NxN) state weighting matrix - 
must use single subscript* 

— 

R 


R 

(UxM) control weighting matrix - 
must use single subscript* 

— 

R 


T 

Interval sampling time 

sec 

R 

. 0833333 


■"To laput a value in the (I,J)th position of the matrix use 'the following algrbraic conversion 
compute appropriate subscript for singly subscripted variables 

K » (J-1) ♦ N + I 

where N is the number of rows of the matrix 


TABLE 2.1-9 
NAM8 INPUT 


VARIABLE 

DESCRIPTION 

UNITS 

TYPE 

IKFAULT 

UMAGl 

Is a 17 x 4 array of steady-state 
values of the systcjm where I 
corresponds to the 1th state and 
J corresponds to the jth run 

- 

R 

68 * 0.0 

UMAG2 

Is a 4 X 4 guidance command array 
where I corresponds to the ith 
state and J corresponds to the 
jth run 


R 

16 * 0.0 

DXZERO 

Is a 17 X 4 array of initial con- 
ditions on the 1th states, con- 
trols and Integrators and to the 
Jth run 


R 

68 * 0.0 
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TABLE 2.1-10 


NAM9 INPUT 


VARIABLE 

DESCRIPTION 

UNITS 

TVPl 

DEFAULT 

TRDT 

Iteration interval for continuous- 
time propagation of Riccatl 
equation 

■ec 

R 

0.5 

SMPLT 

Iteration interval for sampling 
time determination 

■ec 

R 

.001 

UAX 

Maximum number of steps allomd 
for propagation of Riccatl 
equation 

• 

I 

30 

HAXD 

Uaximum number of steps allowed 
for propagation of discrete 
Riccatl equation 


I 

SO 

EPS 

Convergence criterion for both 
discrete and continuous designs 

- 

R 

.01 

IlUX 

Uaximum number of time steps for 
simulation 

- 

I 

300 

IRON 

Number of simulation runs 

- 

I 

1 

BELT 

Iteration Interval for continuous 
simulation 

sso 

R 

.1$ 

ITY1»{1) 

controller type flag 

ITYP(l) “ 1 dynamic compensation 

achieved by control-rate weight 

ITYP(l) > 1 dynamic compensation 

achieved by PII 

ITYP(1) ” 2 velocity control 

system 

ITYP(l) ■ 3 attitude control 
system 


I 

a 

ITYPCa) 

ITYP(2> - 1 UilAGl is used in Sec- 
tion 7.1 

1TYP(2) - 2 UilAGl is not used in 
Section 7.1 

ITYP(2) - 3 punches out nominal 
control and angle position 
1TYP(2) • 4 or 5 determine and 
simulate discrete PI controller 
ITYP(2) “ 4 attitude control 
sy.stem 

1TYP(2) S velocity control 

system 

ITYP{2) “ 6 or 7 determine and 
simulate discrete control-rate 
system 

ITYP(2) - 6 attitude control 
systera 

1TYP(2) “ 7 velocity control 
system 

ITYP(2) - 8 punch out F 6x6 and 
G 6x4 for trim In [AERO 
1TYP(2) “ 9 punch discrete con- 
trol matrix 


I 

1 


* » 


r > 


* * 
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2 . 2 DIGADAPT OUTPUT 

The DIGADAPT output takes two forms: printed reports 

on Fortran logical unit 6 (TAPE6) and punched output on Fortran 
logical unit 7 (TAPE?). Table 2.2-1 contains a summary of 
the printout from the program including a paraphrase, a de- 
scription, the variables and the subroutine associated with 
the report. Table 2.2-2 contains a list of all the variables 
that can be punched for the SCHED regression analysis. It 
includes a description of the variables, the format of the 
punched output, the subroutine generating the output and the 
control flag(s) necessary to turn the punch command "on”. 
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TABLE 2.2-1 


OUTPUT REPORTS FROM DIGADAPT 


REPORT 

DESCRIPTION - VARIABLE 

SUBROUTINE 

CASE NUBIBER = 

Give particular case of aerody- 
namic coefficients chosen - will 
print out only in multiple case 
run starting with second case - 
NCASE 

MAIN 

DIGADAPT 

' Abstract of the program 

MAIN 

INPUT PARAMETERS 

The flags, defaults and/or 
override values that will be 
used in all cases to follow - 
. ICHPTR and other parameters 

MAIN 

WEIGHTING ON STATE 

Values used in computing the 
Q matrix 

MAIN 

WEIGHTING ON INTEGRA- 
TORS FOR VELOCITY- 
COMMAND SYSTEM 

X position, y position, z posi- 
tion, yaw integral - Tl, T2, 

T3, T4 

MAIN 

TRIM VALES 

^B’ "^C’ S’ ^R’ * 

MAIN 

NORMAL EOJ 

Job has terminated as expected 

MAIN 

SYSTEM MATRIX F 

N X N system matrix F 

AERO 

INITIAL VELOCITIES 


AERO 

XDOT, YDOT, ZDOT, THETA 
PHE, PSI 

X, y, z, e, 4>, ip 

FGCOMP 

PO, QO, RO 

Nominal control 

FGCOMP 

NOISE COVARIANCE MATRIX 

QBAR 


STATE VECTOR BOUND 

SVB 

SID 

SAMPLING TIME 

Sampling time T as determined by 
SAMPLE 

SID 

STATE POSITION WHOSE 
BOUND WAS EXCEEDED 

The diagonal value of co- 

variance matrix which was 
greater than state vector bound 

SID 

TIME HISTORY OF DIA- 
GONAL ELEMENTS COVARI- 
ANCE 

Beginning subroutine SAMPLE 
with value of iteration inter- 
val for sampling time - the 
results to follow are U, V, W, P, 
Q,R, THETA, PHI, PSI TIME 

SAMPLE 

STATE WEIGHTING MATRIX 

NxN Q matrix 

GAIN 

CONTROL WEIGHTING MATRIX 

MxM R matrix 

GAIN 

OPTIMAL GAIN MATRIX 

MxN C matrix 

GAIN 
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TABLE 2.2-1 

OUTPUT REPORTS FROM DIGADAPT (Continued) 


REPORT 

DESCRIPTION - VARIABLE 

SUBROUTINE 

RICCATI SOLUTION EQUATION 

NxN X matrix 

GAIN 

DISCRETE RICCATI EQUA- 
TION 

MxN XD matrix 

GAIN 

DISCRETE OPTIMAL GAIN 

MxN CD matrix 

GAIN 

TRANSFORMATION MATRIX 
FOR PI CONTROLLER 

NRxP Q matrix where NR=N-M-P 

GAIN 

FEEDBACK GAIN MATRIX 
FOR PI CONTROLLER 

MxM R matrix 

GAIN 

FEEDBACK GAIN MATRIX 
FOR PI CONTROLLER 

MxNR QHAT matrix 

GAIN 

PROPAGATION OF CONTIN- 
UOUS TIME RICCATI EQUA- 
TION 

TRDT - time step for iteration. 
MAX - maximum of number of iter- 
ations allowed 
EPS - convergence criterion 
NUM - number of iterations 
computed 

RICTCT 

BEGIN SOLUTION OF ALGE- 
BRATIC RICCATI EQUATION 

Use Kleinman scheme with con- 
vergence criterion 10~5 

TMRIC 

ALGEBRAIC RICCATI 
SOLUTION IN 

Give number of iterations to 
satisfy convergence criterion 

TMRIC 

RICCATI SOLUTION IS PSD 

Riccati solution is positive 
semi-definite 

TMRIC 

RICCATI NOT CONVERGED 

Computation continues as if con- 
vergence occurred - value of 
solution matrix x is output 

TMRIC 

RICCATI BLOWUP 

Solution X has diverged 

TMRIC 

LINEAR EQN ALGORITHM 
NON-CONVERGENT 

Gives number of Iterations, IT, 
attempted to solve linearized 
Riccati equation 

MLINEQ 

NO RICCATI SOLUTION 

Gives number of iterations, I, 
attempted to solve Riccati 
equation 

TDREQ 

PROPAGATION OF DISCRETE - 
TIME RICCATI EQUATION 

MAXD - maximum number of itera- 
tions allowed 

EPS - convergence criterion 
I - number of iterations 

TDREQ 

BEGIN SOLUTION OF 
STEADY-STATE DISCRETE 
MATRIX RICCATI EQUATION 

Use Kleinman scheme with con- 
vergence criterion of 10"^ 

DRIC 
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TABLE 2.2-1 

OUTPUT REPORTS FROM DIGADAPT (Continued) 


REPORT 

DESCRIPTION - VARIABLE 

SUBROUTINE 

STEADY- STATE DISCRETE 
RICCATI SOLN 

Gives number of iteration, IT, 
performed 

DRIC 

RICCATI SOLUTION IS PSD 

Riccati solution is positive 
semi-definite 

DRIC 

RICCATI NON-CONVERGENT 

Computation continues despite a 
non-convergent solution 

DRIC 

RICCATI BLOWUP 

Solution diverging 

DRIC 

STEADY- STATE MATRIX FOR 
GUIDANCE SYSTEM S12 

DUMl 

PID 

STEADY CONTROL MATRIX 
FOR A CONTROL SYSTEM S22 

DUM3 

PID 

COVARIANCE MATRIX FOR 
PROCESS NOISE 

SIGMAD 

EOD 

COVARIANCE MATRIX FOR 
OBSERVATION NOISE 

THETAD 

EOD 

OBSERVATION MATRIX FOR 
DISCRETE SYSTEM 

MEAT 

EOD 

RICCATI EQUATION STEADY 
STATE SOLUTION 

XD 

EOD 

KALMAN GAIN FILTER 
MATRIX 

DT 

EOD 

DISCRETE CONTROL MATRIX 

XMLIH 

EOD 

EIGENVALUES OF ACL 

Closed- loop system 

EIGEN 

EIGENVALUES OF FCL 

Closed- loop discrete Kalman 
filter system 

EIGEN 

EIGENVALUES 

Printout eigenvalues, period, 
damping ratio, sunplitude and 
frequency 

FREE 

SYMBOL MEANING AND UNITS 
OF SIMULATION 

Gives name of symbols used in 
time history calculation 

TIME 

STEADY- STATE VALUES 

UMAGl 

TIME 

GUIDANCE COMMAND MATRIX 

UMAG2 

TIME 

CONTINUOUS- TIME SIMULA- 
TION PII 

If ITYP(l) equals 2 

TIMHIS 

CONTINUOUS-TIME SIMULA- 
TION OF DYNAMIC CON- 
TROLLER 

If ITYP(l) equals 1 

TIMHIS 
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TABLE 2.2-1 

OUTPUT REPORTS FROM DIGADAPT (Continued) 


REPORT 

DESCRIPTION - VARIABLE 

SUBROUTINE 

WRITE OUT TIBffi HISTORIES 
U,V,W,P,Q,R, THETA, PHI 
PSI,DELTA-B, DELTA-S, 
DELTA-R,V1,V2,V3,V4, 
DXDOT , DYDOT , DZDOT , U1 , U2 , 
U3, UDOT,VDOT,WDOT, 
DELBDOT, DELCOT.DELSDOT, 
TIME 

PX matrix is written after each 
run (maximum of 4) is completed 

TIMHIS 

STEADY STATE HAS NOT 
BEEN ATTAINED 

Gives ith variable of PX for 
which absolute difference is 
greater than 10”'*. 

STEADY 

STEADY STATE 

Writes out last time step value 
for PX matrix, SS, on the L^“ 
component . 

STEADY 

PEAK VALUE 

Writes out maximum value for i^^ 
component of PX=PV 

STEADY 

TIME TO REACH NINETY 
PERCENT OF STEADY STATE 

Time to reach 90% of steady 
state - T90 

STEADY 

TIME TO REACH A FIVE 
PERCENT BOUND OF STEADY 
STATE 

T05 

STEADY 

DISCRETE SIMULATION OF 
PI I INTEGRATOR 

Begin subroutine DTMHIS 

DTMHIS 

WRITE OUT TIME HISTORIES 
(SAME AS TIMHIS) 

DISCRETE SIMULATION OF 
PI 

PX is written after each run 

DTMHIS 

ITYP (2) equals 4 or 5 

PIDHIS 

DISCRETE SIMULATION 
OF DYNAMIC CONTROLLER 

ITYP (2) equals 6 or 7 

PIDHIS 

WRITE OUT TIME HISTORIES 

PX is written after each run 

PIDHIS 

Holerith title 

In print matrix routine hole- 
rith title is printed along 
with number of rows and columns 

PRMR 

FATAL ERROR - SINGULAR- 
ITY DISCOVERED 

The propagation of covariance 
matrix has become singular - 
W1 , 

PROP 

SINGULARITY DISCOVERED 

Singularity has been found in 
Kalman filter update - W2 

SMPL 

FATAL ERROR - MEASUREMENT 
NOISE MATRIX IS SINGULAR 

R 

RSET 
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TABLE 2.2-2 

PUNCHED OUTPUT FROM DIGADAPT 


VARIABLE 

DESCRIPTION 

FORMAT 

SUBROUTINE 

CONTROLLING* 

FLAG 

F 

Upper 6x6 matrix of aircraft 
system matrix used in dynamic 
trim computation and the upper 
4x6 pseudo-inverse of aircraft 
control system matrix used in 
dynamic trim computation 

6E12.5 

AERO 

ITyP(2)-8 

NCASE 

Case number corresppndlng 
to aerodynamic coefficients 

15 

FCiCOMP 

NPUNCH-1 

XDOT.YDOT 
2D0T, THETA 
PHE.PSI 

X. y. z, 0, iJ* 

6E12.5 

FGOOMP 

ICHPTR(2,1)-1 

PS I DOT 

♦ 

E12.5 

FGCOMP 

ICHPTR(2,1)-1 

PO,QO,RO 

Nominal control 

3E12.S 

FGCOMP 

ITYP(2)«3 

ICHPTR(2,1)-1 

UO,VO,WO 
THETA, PHE, 
PSI 

“O' '^0' '*0* 

6E12 . 5 

PGCOMP 

ICHPTR(2,1)-1 

c 

mxn optimal gain matrix 

6E12.5 

GAIN 

ICHPTR(4,3)-1 

CD 

mxn discrete optimal gain 
matrix 

6E12.5 

GAIN 

ICHPTR(4,5)-1 

R 

Gain matrix for delayed error 
state 

6E12.5 

GAIN 

ICHPTR(4,6)=1 

RHAT 

Gain matrix for current error 
state 

6E12.5 

GAIN 

iaiPTR(4,6)-l 

QHAT 

State feedback gain matrix 

6E12.5 

GAIN 

ICHPTR(4,6)-1 

DT 

(9xL) Kalman filter gain matrix 

6E12.5 

EOD 

ICHPTR(5,1)-1 

FWORK 

Discrete 9x9 discrete system 
matrix 

6E12 . 5 

EOD 

ICHPTR(5,1)-1 

ITYP(2)=9 

XMLIH 

Discrete 9xra control matrix 

6E12 . 5 

EOD 

ICHPTR(5,1)-1 

ITYP(2)-9 


♦In addition to the controlling flags listed HPUNCH must be set to one in the header input card 
for punched output. 
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3 . DETAIIiED DESCRIPTION OF DIGADAPT 


DIGADAPT is comprised of seven program chapters : the 
main executive, system matrix computation, sampling interval 
determination, linear-optimal gain design, calculation of 
eigenvalues, and time history simulations. This chapter dis- 
cusses the significant subroutines of DIGADAPT and provides 
flowcharts and tables of input and output. Figures 3-1 to 3-3 
gives a structural overview of the relationships between rou- 
tines in outline form. A brief discussion of the service 
utilities also is included. 


3.1 MAIN EXECUTIVE 

The main program reads the data for a given run and 
calls all subsequent program chapters. On the first case for 
a given execution of DIGADAPT, the program reads a header card 
and then reads data via NAMELIST to override the BLOCK DATA 
initialization. On subsequent cases of a multiple case run, 
only the header card is needed (see Fig. 2.1-1). The program 
terminates when the last case input on Fortran logical unit 5 
has been processed. 

The Output from the main executive consists of an 
abstract of the DIGADAPT program along with a listing of ini- 
tialization parameters (printed out for the first case only). 
Figure 3.1-1 gives an overview of the main executive. 
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Overview of Program DIGADAPT 
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Subroutine TIME and Service Routines 
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Figure 3.1-1 Flowchart of Main Executive (Segment 2 of 2) 


3.2 SYSTEM MATRIX COMPUTATION 

The purpose of subroutines AERO and FGCOMP is to place 
the stored aerodynamic coefficients into the correct positions 
of the system matrix F. The ordering of the states for this 
program chapter is as shown in Table 3.2-1. 
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TABLE 3.2-1 
ORDER OF STATES 


STATE 

VARIABLE 

DESCRIPTION 

u 

BODY-AXIS VELOCITY - POSITIVE FORWARD 

V 

BODY-AXIS VELOCITY - POSITIVE RIGHT 

w 

BODY-AXIS VELOCITY - POSITIVE DOWN 

e 

INERTIAL-BODY PITCH EULER ANGLE 


INERTIAL-BODY ROLL EULER ANGLE 


INERTIAL-BODY YAW EULER ANGLE 

P 

BODY- AXIS ROLL RATE 

q 

BODY-AXIS PITCH RATE 

r 

BODY-AXIS YAW RATE 


DIFFERENTIAL COLLECTIVE 


GANG COLLECTIVE 


GANG CYCLIC 

<5r 

DIFFERENTIAL CYCLIC 


Appendix A describes how the disk storage for the 
aerodynamic coefficients is prepared and describes how to 
choose appropriate calls. This disk file must be created 
before running DIGADAPT. Figure 3.2-1 gives the flowchart of 
AERO, and Table 3.2-2 gives the inputs and outputs for FGCOMP. 


3.3 SAMPLING INTERVAL DETERMINATION 

The purpose of subroutines SID and SAMPLE is to aid 
in the determination of the sampling time of a continuous, 
linear-time-invariant system. In this program chapter, the 
system matrix F has three added gust states. The sampling 
time is found by propagating the covariance matrix equation 

X(t) = FX(t) + X(t)F'^ + Q 
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Figure 3.2-1 Flowchart of Subroutine AERO 


TABLE 3.2-2 

INTERFACE TO ROUTINE FGCOMP 

INPUTS 

N = number of states 
M = number of controls 
NCASE = case number 

NPUNCH = 0 if no punch output is desired 
1 if punch output is desired 

OUTPUTS 

F = NxN state matrix 
G = NxM control matrix 
UO = body velocity along body x-axis 
VO = body velocity along body y-axis 
WO = body velocity along body z-axis 
XDOT - body velocity along inertial x-axis 
YDOT = body velocity along inertial y-axis 
ZDOT = body velocity along inertial z-axis 
ANGLE = 3x3 transformation matrix from inertial 
coordinates to body coordinates 
ANGLI =3x3 transformation matrix from body to 
inertia} coordinates 
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where F is the augmented system matrix, X is the state covari- 
ance* Q is the covariance matrix of the system excitation noise, 
and X(0) = 0. The sampling time is determined when any diago- 
nal element of X exceeds the input state vector bounds (SVB). 
Figure 3.3-1 shows how this is accomplished and Table 3.3-1 
gives the interface of SAMPLE with SID. 

Rat7t» 



Figure 3.3-1 Flowchart of Subroutine SID 
*X is equivalent to engineering notation P. 
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TABLE 3.3-1 

INTERFACE TO ROUTINE SAMPLE 


INPUTS 

N = dimension of system matrix 
SVB = state; N-vector bound 
QBAR = NxN noise covariance matrix 
F = NxN system matrix 

OUTPUTS 

T = sampling time 

K = the number of the state whose state 
vector bound was exceeded 


3.4 LINEAR-OPTIMAL GAINS 

The purpose of this program chapter is to design 
various linear-optimal gains (Fig. 3.4-1). GAIN is divided 
into continuous and discrete modes of design ^ with and with- 
out integrators. 

The following equations show the controllers which 
can be simulated, where the K's shown are continuous gains, 
the C's are discrete gains, and $ and F are discrete equiva- 
lents of F and G. 

Continuous Time Dynamic Controller 

x(t) = Fx(t) + Gu(t) 

u(t) = KgUCt) + K^x(t) + Ly^ 

Continuous Time PII Controller 

x(t) = Fx(t) + Gu(t) 

u(t) = KgUCt) + K^xC-t) + KgZCt) 

i(t) = Tx(t) - y^ 
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Figure 3.4-1 Flowchart of Subroutine GAIN 

(Segment 1 of 3) 
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Figure 3.4-1 Flowchart of Subroutine GAIN 

(Segment 2 of 3) 


REPRODUCIBILnY 
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Figure 3.4-1 Flowchart of Subroutine GAIN 

(Segment 3 of 3) 


Discrete PI I Controller 
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Simplified Discrete PII Controller 


~x" 

u 

S5 

■ $ r 0" 
0 10 


'x! ;3| 

1 

+ 

0 

Atl 

[-Cl -s -s] 

”x~ 

u 

_Z. 

k+1 

_AtT 0 I_ 


UJ 

k 

0 


.Z. 


0 

0 

At I 


id 


Simplified Discrete Dynamic Controller 
Sk+l = 2k 

ak+1 “ Sk - ’’^k> + 

Simplified Discrete PI Controller 
Sk+l = 2k * '■^k 

% = Hk-1 - ’fSk.i) + - TXk) + Cg(x^) 


3.4,1 Q and R Computation 


The purpose of QRCXDMP is to compute the state weight- 
ing matrix Q and the control weighting matrix R. Note that if 
NSIGMA (input variable of header card) is set to one then the 
value of SIGMA(I), I = 1 to 16 will default to the BLOCK DATA 
value; alternatively, NAMELIST may be used to override the 
default value. Note that the BLOCK DATA values default to 
negative values; hence, the program sets SIGMA(I), 1=1 to 
16 to zero. The SIGMA values are then used to compute the 
diagonal terms of the Q matrix. The diagonal values of the 
R matrix are computed by taking the inverse squared of the 
appropriated maximum control rate. Table 3.4-1 gives the 
interface between QRCOMP and GAIN. 
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TABLE 3.4-1 

INTERFACE TO ROUTINE QRCOMP 


INPUTS 

N B Dumber of states 
11 number of controls 
NSIGUA * 0 If SIGMAS are to be calculated 
1 if SIGMAS are to be input 
NDOT > 0 If diagonal values of the R matrix 
are to be computed 
1 if diagonal values of the R matrix 
are to be input 

OUTPUTS 

Q • NxN state weighting matrix 
R ■ HxM control weighting matrix 


3.4.2 Guidance Commands 


The purpose of TCOMP is to augment the states of the 
system matrix, F, the control matrix, G, and the state weighting 
matrix Q. The added states have the form ^ = Tx - UMAG2 where 
UMAG2 are the guidance commands, and T is the transformation 
matrix converting x to UMAG2 states. Table 3.4-2 shows the 
interface between TCOMP and GAIN. 


TABLE 3.4-2 

INTERFACE TO ROUTINE TCOMP 


INPUTS 

N * number of states without T 
H > number of controls 
P “ NxN state matrix 
F “ NxU control matrix 
Q * NxN state weighting matrix 
ITYP - type of guidance used 
Tl) 

weighting on extra states for attitude system 

V*} 

OUTPUTS 

N • number of states with T 

P ■ number of states in extra, states added to system 
F ■ NxN new state matrix 
G “ NxN new control matrix 
Q «■ NxN new state weighting matrix 
angle > 3x3 transformation matrix from Inertial 
coordinates to body coordinates 
ANGLI * 3x3 transformation matrix from body to 
inertial coordinates 
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3.4.3 Continuous Regulator Gains 

REG designs an optimal linear regulator by finding 
the steady-state solution to the Riccati equation 

p = PF + f’^P + q - PGR“^ G*^P 

In addition, the optimal feedback gain K and the closed-loop 
system matrix ACL are calculated. REG first propagates the 
Riccati equation (using subroutine RICTCT) from P = 0 until 
either the maximum number of iterations, MAX, is reached or 
the convergence criterion, EPS, is satisfied. It then takes 
the last value P from RICTCT as an Initial input guess for 
TMRICE, an entry point in TMRIC. Table 3.4-3 shows the inter- 
face between REG and GAIN. 

TABLE 3.4-3 

INTERFACE TO ROUTINE REG 

INPUTS 

N = number of states 
M = number of inputs 
A = NxN system matrix F 
B = NxM input matrix G 
Q = NxN state weighting matrix 
R = MxM control weighting matrix 

OUTPUTS 

X = NxN matrix Riccati equation solution 
G = MxN optimal gain matrix 
ACL = NxN optimal closed-loop matrix 


The purpose of RICTCT is to solve the Riccati equation 
-P(t) = P(t)F + F^^P(t) + Q - P(t) GR"^GP(t) 
using backward integration (Ref , 1) . Consider the equations 
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If a solution exists on the interval (t,t^), then the solu- 
tion has the property that 


P(t) = Y(t) X'^(t) 


For computer implementation, the system matrix is discretized 
by 




0 ( t , T ) = exp [ 

-Q -F 

where 0 (t,x) is partitioned into four nxn submatrices 


F -GR"^g'^^ 


(t-T) 


0 = 


11 

®12 

21 

®22 


If t-T is set to a constant value At, then the solution for 
P can proceed as follows: 

P(tjj) = Pq 

p((n-l)At) = [ 02i(At) + 022(At)P(nAt)] [ ©^^^(At) 

+ P(nAt)] 

This computation is implemented by calls to RSET which compute 
the transition matrices and RCT v/hich solves the Ricatti equation. 

Finally, upon leaving RCT, if NUM=MAX, no steady- 
state solution has been reached; if NUM<MAX a steady-state 
solution has been reached, where NUM is the number of iterations 
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performed. In addition the values of EPS, MAX, NUM and TRDT 
are printed. Then the closed-loop system ACL matrix is com- 
puted by r - GR~^g'^P. Table 3.4-4 gives the interface be- 
tween RICTCT and REG. 


TABLE 3.4-4 

INTERFACE TO ROUTINE RICTCT 


INPUTS 

F = NxN system matrix 
G = NxM control matrix 
Q = NxN state weighting matrix 
R = MxM control weighting matrix 
TRDT = time step of propagation 
MAX = maximum nvimber of time steps allowed 
EPS = criterion for reaching steady state 

OUTPUTS 

X = NxN Riccati solution at time NUM*TRDT 
C = X = MxN optional feedback gain 

ACL = F - G*C = NxN closed- loop system 


The purpose of TMRIC is to solve the algebraic Riccati 


equation 


0 = PF + F^P + Q - PGR"^GP 


using Newton's Method. Let Pj^, K=0,1,... be the (unique) posi- 
tive semidefinite solution of the linear equation 


0 = + Q + lTrl^ 

where, recursively, 

Lj^ = R~^g'^Pj^_^ k=l,2,... 

F, = F - GL, 
k k 

Lq is chosen so that the matrix F - GL^ has eigenva,lues with 
negative real parts. Thus 

0 < P < Pk+i-^^k - • • -^0 k=U,l,2, ... 

and 
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The above scheme linearizes the Riccati equation and can be 
shown to have quadratic convergence. 


The linearized form of the Riccati equation is solved 
by subroutine MLINEQ with the following decisions made in 
TMRIC 

• Has the solution converged? 

• Has the solution diverged? 


• Has the 
Table 3.4-5 shows the 


solution not diverged and not converged? 
interface between TMRIC and REG. 


TABLE 3.4-5 

INTERFACE TO ROUTINE TMRIC 


INPUTS 

N = dimension of square arrays 
A = NxN system matrix 
S = NxN matrix = G*R~1 *gT 
Q = NxN state weighting matrix 

OUTPUTS 

X = NxN positive-definite solution of matrix 
Riccati equation 

Z = A-S*X = NxN closed-loop matrix 


3.4.4 Calculation of PHI and GAMMA 


The purpose of QRMHAT is to discretize a linear-time- 
invariant continuous system and its associated cost function 

X = Fx + Gu 


J(u) = 


T T 

X Qx + u Ru dt 


In its discretized form, the first equation becomes 


2k+l “ * ’’Hk 
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Similarly, the second equation becomes 


J(u) = I + 2x^ Mu^ -k 

k“0 

Upon discretization, subroutine QRMHAT produces the following 
results (computer sjmibols are used on the left, engineering 
notation on the right): 


PHI = I 
k= 

rT “ 


(FT)^ _ .FT 
k ' ® 

k=0 


MHAT = 


QHAT 


F-s J 
e * Q{ 

L' o Vo 


^ ■ j; j. ■ 
■i: 


fT 


o 


Ft , 
e dT 


F*^T Ft 
e^ ^Qe dT 


e^''^ dwjds 


RHAT = TR + G*^ 

f1 

r® f\ 

e dw 

Q • 

e^™dw 


o( 

-■'o 


. O 


ds )G 


This computation is completed by EXPI computing PHI, by SETA 
computing QHAT, and by ALPHA computing RHAT and MHAT, Table 
3.4-6 shows the interface between this program and GAIN. 

The purpose of MQCOST is to find equivalent matrices 
FEQ and QEQ so that the optimal regulator problem 


and 


J = 


X = Fx + Gu 


(x*^ Qx + 2x*^ Su + u"^ Ru) dt 


is similar to 
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TABLE 3.4-6 

INTERFACE TO ROUTINE QRMHAT 


INPUTS 

M = number of controls 
N = number of states 
T = sampling time 
F = NxN state matrix 
G = NxM matrix 
Q = NxM matrix 
R = MxM weighting matrix 

OUTPUTS 

PHI = NxN discretized system matrix 
DT = NxM discretized output matrix 
QHAT = NxM discretized state cost matrix 
MHAT = NxM discretized cross coupling cost matrix 
RHAT = MxM discretized control cost matrix 


X = (FEQ)x + Gu 

and 

fCO 

J = I (x"^(QEQ)x + u*^Ru)dt 
' o 

Table 3.4-7 shows the interface between MQCOST and TDREQ. 


TABLE 3.4-7 

INTERFACE TO ROUTINE MQCOST 


INPUTS 

N = number of states 
M =‘ number of controls 
F = NxN matrix 
Q = NxN matrix 

OUTPUTS 

FEQ = NxN equivalent F matrix = F - G*R“l*sT 
QEQ = NxN equivalent Q matrix = Q - S*R-l*sT 
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3.4.5 Discrete Regulator Gains 

The purpose oi DREQ is to design an optimal linear 
regulator for the discrete time, linear system 


^k+i = '■^k 


according to the quadratic cost function 


__ 

^ ' klo 


The regulator is defined by the optimal control Uj^ = -Kxj^ where 
and 


K = (R + r'^pr)*^ r’’’p4i 


^k+l “ (^ci')ik 

For DREQ the cost function is 


and 


ou 

J ^ ^-k ®-k -k^-k^ 

k“ 0 

T -1 T T 

K = (R + r^pr) (r p$ + s ) 


Subroutine DREQ insures that the first P used in DRIC 
is chosen correctly. The routine propagates the Riccati solu- 
tion P (using subroutine TDREQ) from P = 0 until either EPS is 
satisfied or MAXD is reached. The solution P, an output from 
TDREQ, is then used by subroutine DRIC, where a final solution 
is obtained.* Table 3 . 4-8 shows the interface between GAIN and DREQ, 


Subroutine TDREQ propagates the discrete-time Riccati 


equation 


\ Pk+i» ^ Pk+i 


♦Note the program variables A, B, X, and G correspond to the 
above variables $,T,P,K respectively. 
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TABLE 3.4-8 

INTERFACE TO ROUTINE DREQ 


INPUTS 

N = number ol states 
M = number of Inputs 
A = NxN system matrix 
B = NxM input matrix 
R = MxM control weighting matrix 
Q = NxN state weighting matrix 

S = NxM control and state crossed weighting matrix 
OUTPUTS 

X = NxN matrix Riccati equation solution 
G = MxN optimal gain 
ACL = NxN optimal closed-loop matrix 


where the variables, PHI, XD, QHAT, RHAT, DT are equivalent 
to variables $, P, Q, R, r respectively. 


This calculation is accomplished by looping on the 
routine SMPL and PRPA, which updates and propagates the matrix 
P. When the convergence criterion EPS is satisfied, as when 

MAXD, the maximum number of iterations is reached, the closed- 

loop system ACLD and the feedback gain CD are computed using 

the equations: 

CD = (R + rpr^^)"^ (r'^p$) 


ACLD = $ - r(CD) 

Table 3.4-9 gives the interface between TDREQ and DREQ. 


The purpose of subroutine DRIC is to solve for the 
steady-state solution of the discrete-time Riccati equation 
by Newton's method. This is achieved by the following algorithm 
Pj^(k=0,l,2, . . . ) is the solution of the linear equation 
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TABLE 3.4-9 


INTERFACE TO ROUTINE TDREQ 


INPUTS 

N = number of states 
M = number of controls 
PHI = NxN discrete system matrix 
DT = NxM discrete output matrix 
RHAT = MxM discrete weighting on controls 
QHAT = NxN discrete weighting on states 
MHAT = NxM discrete cross coupling weighting 

OUTPUTS 

XD = NxN final propagation solution to the 
Riccati iteration 
CD = MxN discrete feedback gain = 

(RHAT + DTT*XD*DT)-^*(DTT*XD*PHI) for TDREQ 
or (RHAT + DTT*XD*DT)-1*(DTT*XD*PHI+MHATT) 
for TDREGl 

ACLD = NxN closed -loop system = PHI - DT+CD 



$3 P, + Z, 
k k k k 


where, recursively 


—1 T -1 

and where P^ is chosen so that the matrix (1 + TR r Pq) $ 
has eigenvalues inside the unit circle. Table 3.4-10 shows 
the interface between DRIC and DREQ. 

3.4.6 Proportional- Integral Controller 

PID computes the gains for the proportional-integral 
controller using the following set of equations. If ITYP(l) = 
4 or 5, then 


(I + TR k=l,2, . . 


»k Pk-i \ ^ « 
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TABLE 3.4-10 


INTERFACE TO ROUTINE DRIC 


INPUTS 

N = dimension of system 
A = NxN system matrix - 
S = NxN matrix = B*R-1*B 
Q = NxN weighting matrix 

X = NxN initial guess matrix for calling DRICl 
TR = NxN gain matrix 

OUTPUTS 

X = NxN positive definite solution of matrix 
Riccati equation 

Z = NxN closed-loop system matrix = (I+S*X) *A 
IOC = -1 if anything goes wrong 


Ft 

--k+1 ~ ® ^ ^ 


e^® ds 1 Gu 


(3.6-1) 


-k " -k-1 r((UMAG 2)^..- Qx^-i) + (RHAT)((UMAG2) - 

+ (QHAT)x^ (3.6-2) 

If ITYP(l) = 6 or 7, then Eq. (3.6-1) remains the same, 
but Eq. (3.6-2) becomes 


-k " -k-1 ^ ^k-1 (RHAT)(UMAG2 - QXj^_^) + (QHAT) 


(3.6-3) 


Table 3.4-11 gives the interface between PID and GAIN. 
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TABLE 3.4-11 
INTERFACE TO ROUTINE PID 


INPUTS 

N = dimension of composite system 
M = dimension of controller 
P = dimension of guidance system 
KPUNCH = logical unit number for punched output 
NPUNCH = 0 if no punched output is desired 
T = sampling time of discrete system 
NR = dimension of system = N-M 
F = NxN system matrix 

OUTPUTS 

R = MxP gain matrix for delayed error state 
RHAT = MxP gain matrix for the current error state 
QHAT = MxNR state feedback gain matrix 

Q = I^xNR trSusf ormation matrix of control system . 


3.5 . ESTIMATOR - OBSERVER DESIGN 


The purpose of this program chapter, subroutine EOD, 
is to perform the estimator observer design calculations 
(Fig, 3.5-1). 

The routine ESTOBS is called to find the steady-state 
discrete-time Kalman filter gains for a continuous -time sys- 
tem. The system matrix F is discretized with interval T. 

Table 3.5-1 shows its interface with EOD. 

Finally ESTOBS calls KFD, whose purpose is to design 
a steady-state discrete Kalman filter. Table 3.5-2 shows the 
relationship between ESTOBS and KFD. 
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No 


Calculit* 
obieivation 
tiMtrix for disa«tc 
iysttm uiing KFLAG 
input to determint 
number of observers 



Covariance 
matrix SIGMAO 
for process noise 



Covariance matrix 
THETAD for 
observation noise 



ESTOBS 


Estimator 
Observer 
corn potation 


Discrete system 
matrix e^' 



RIccatl eq. steady 
state solution XD 



Figure 3.5-1 Flowchart oi Subroutine EOD 

(Segment 1 of 2) 
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TABLE 3.5-1 

INTERFACE TO ROUTINE ESTOBS 


INPUTS 

N = dimension of system 
L = dimension of observation 
F = NxN system matrix 
H = NxN observation matrix 
SIGMAD = NxN process noise covariance matrix 
THETAD = LxL observation noise covariance matrix 
(dimensioned at least NxN) 

EPS = convergence criterion for TDREQ 
MAXD = maximum number of iterations for TDREQ 
T = sampling time for filter 

OUTPUTS 

S = 9x9 discrete time steady-state Riccati 
equation solution 

DK = 9xL discrete time Kalman filter gain 
FCL = 9x9 discrete Kalman filter closed-loop 
matrix = F-K*H 


TABLE 3.5-2 

INTERFACE TO ROUTINE KFD 


INPUTS 

N = niimber of states 
L = number of observed variables 
A = NxN system matrix 
H = LxN observation matrix 

R = LxL observation noise covariance matrix 
Q = NxN process noise covariance matrix 
D = NxL correlated noise matrix 

OUTPUTS 

S = NxN discrete Riccati equation steady 
state solution 

K = NxL Kalman filter discrete gain 
FCL = A-K*H 
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3.6 EIGENVALUE COMPUTATION 

The purpose of this program chapter, subroutine EIGEN 
is to calculate the eigenvalue of a given matrix. Consider Fig. 
3.6-1, which shows that the eigenvalues of three matrices are 
of interest to DIGADAPT. They are ACL, the closed loop system, 
ACLD, the closed-loop discrete system, and FCL, the closed-loop 
Kalman filter system. Note that if the user so desires, he can 
calculate the eigenvalue of any matrix by a direct call to FREE, 
where it should be noted that, in calculating the eigenvalues, 
the matrix of interest is destroyed. 

The purpose of FREE is to calculate the eigenvalues 
of various modes of motion of the given aircraft system matrix, 
along with the corresponding natural frequencies, damping 
rates, periods, and times to half amplitude. In turn, FREE calls 
two programs HSBG and ATEIG, which have been incorporated from 
the System/360 Scientific Subroutine Package (Ref. 4). 

The subroutine HSBG reduces an NxN real matrix A 
by similarity transformation to upper almost-triangular 
(Hessenberg) form. Each row is reduced in turn, starting from 
the last one by applying a suitable right elimination matrix, 
and similarity is achieved by also applying the left inverse 
transformation. Thus the eigenvalues are preserved. Follow- 
ing the subroutine HSBG, ATEIG computes the eigenvalues of a 
real upper almost-triangular matrix using the double-QR 
iteration of J.G.F. Francis (Ref. 4). Table 3.6-1 shows the 
interface between EIGEN and FREE. 

3.7 TIME HISTORY SIMULATIONS 

The purpose of this program chapter, subroutine TIME, 
is to control time history calculations for continuous-time 


R^lPKODUClBlLrrY OF THE 

i T -mi. m;' TS PHOT?. 
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Figure 3.6-1 Flowchart of Subroutine EIGEN 

(Segment 2 of 2) 


TABLE 3.6-1 

INTERFACE TO ROUTINE FREE 


INPUTS 

N = number of state in the system 
F = NxN aircraft system matrix of NxN closed- 
loop system matrix 

OUTPUTS 

RR = n-vector with real parts of the 
eigenvalues of F 

RI = n-vector with imaginary parts of the 
eigenvalues 
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3.7.1 Linear Time Invariant Calculation 

TIMHIS simulates a continuous linear time invariant 
(LTI) system 


X 

u 

ij 


= [acl] 


V 


"0" 


'O' 

u 

- 

0 

rUMAG2l + 

I 

X 


.1. 

L. J 

-0. 


[uMAGlj 


using state-transition-matrix propagation of the state. The 
propagation equation is 


-k+1 


^k 

^k+1 

= PHI 

^k 

^k+1 


^k 


+ FgCDELU) + r^(UMAGl) 


where 


PHI = e 


(ACL) t 




^2 


/•At 


• At 


e^ACL)s 


e^ACL)s 


0 

I 

0 

0 

0 

-I 


For each run (up to a maximum of 4) there is a maxi- 
mum of IMAX time steps; for each time step 31 variables are 

■ « 

il, V, *, 6 g, 

Table 3.7-1 shows the interface between TIMHIS and TIME. 


calculated including u, w, p, q, r, 6, Tp , c|) , 6g, 6^, 


^1’ ^2’ ^3' ^4' '^x ' ' 


U2/ Ug, u^, 


S’ 

G’ ^R' 
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TABLE 3,7-1 


INTERFACE TO ROUTINE TIMHIS 


INPUTS 


N 

wz 

total niunber of states and controls 

U 

m 

number of controls 

P 


number of commands 

NR 


N - M - P 

ACL 

X= 

NxN closed loop matrix 

G 

= 

NxM control matrix 

UO) 

VO 

= 

Initial body velocity 

W0» 

I RUN 

= 

number of time histories up to 4 

UUAG2 

s=. 

Px4 guidance inputs to the system 

UHAGl 

= 

Nx5 feedforward control input - 

C 

x: 

HxN feedback gain of system 

DELT 

«= 

interval of discretization 

ITYP(l) 

B 

1 no integrator states 

ITYP(l) 

= 

2 velocity guidance system 

ITYP(l) 

= 

3 attitude guidance system 

ITYP(2) 

= 

2 UMAGl is not used and DELU is. set to zero 

ITYP(2) 

= 

1 UMAGl is turned on 

DELU 

= 

C*UMAG1 

DELV 

= 

UMAG2 

lUAX 

= 

maximum number of iteration steps 

DXZERO 

mm 

Nx5 initial condition of the states, con- 
trols and integrators 

OUTPUTS 

PX 

ss 

see discussion of output reports in Table 
2.2-1 


3.7.2 Discrete PII Calculation 


The purpose of DTMHIS is to simulate a discrete-time 
system. The system simulated is 



The exact system is simulated if ICHPTR(7 j 3)=1; an approxi 
mation to the discrete-time system is simulated if 
ICHPTR(7 , 4)=1 . The approximate equation is 
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Table 


3 

the PI 
lation 


u 






k+1 


T>HI 

0 

tt 


DT 

I 

0 


IJU 



■ 0 ■ 
0 

TI 


[uMAG2]j 


This routine uses the same run logic as TIMHIS. 
3.7-2 shows the relationship between DTHMIS and TIME 


TABLE 3.7-2 

INTERFACE TO ROUTINE DTMHIS 


INPUTS 

N = number of states 

If = number of controls 

NR = number of aircraft states 

P = number of integrator states 

T * sampling time 
F = NxN system matrix 
d = MxN feedback gains matrix 
UMAG2 = Pxl set point or guidance command 
DXZERO = Nxl initial' condition of states 
ACLD » control discrete input matrix 

OUTPUTS 

PX = see discussion of output reports 
in Table 2.2-1 


7.3 Discrete PI and Control Rate Calculation 


The purpose of PIDHIS is to perform a simulation for 
control system. If ITYP(2) = 4 or 5, then the simu- 
uses 

-k+1 

= + R(UMAG2 


$Xk + ruj^ 


” *^-k-l^ RHAT(UMAG2 - Qxj^) + QHAT x^ 
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where PHI 


■c :) ^ 


FT 

$ = e and r = 


fT 


o 


e''® dsl , G 


If ITYP(2) = 6 or 7, then the program simulates 


^k+l = » 5k Hk 


^ik ' ^k-1 


+ RHAT(UMAG2 - + (QHAT)Xj^_j^ 


PIDHIS also has the same time history simulation 
logic as TIMHIS and DTMHIS. Table 3.7-3 shows the interface 
between PIDHIS and TIME. 


TABLE 3.7-3 

INTERFACE TO ROUTINE PIDHIS 


INPUTS 


PHI 

r: 

NxN discrete 

R 

B 

MxP matrix 

RHAT 

= 

yixP matrix 

DXZERO 


Nx4 = initial condition of states 
and controls 

UUAG2 

= 

Pxl guidance command 

QHAT 

= 

MxNR matrix 

Q 

= 

HxNR command transformation matrix 

N 

B 

dimension of system + control 

P 

B 

dimension of guidance command 

NR 

= 

N-M 

T 

B 

sampling time 

1 OUTPUTS 

PX 

ss 

see discussion of output reports in I 

1 Table 2.2-1 J 


3.8 SERVICE UTILITIES 

A number of utility routines are used by the sub- 
routines described in earlier sections. Table 3,8-1 lists 
these service programs, gives an appropriate description of 
each, and make references to Tables 3.8-2 through 3.8-33 that 
describe their interface. 
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TABLE 3.8-1 
UTILITY ROUTINES 


SUBROUTINE 

DESCRIPTION 

TABLE 

OF 1/0 

DOT 

Take dot product between two vectors 

Table 

3. 8-2 

CONTRL 

Calculate control and guidance 
inputs 

Table 

3.8-3 

STEADY 

Calculate steady state in a time 
history simulation 

Table 

3.8-4 

HSBG 

Transform matrix into Hessenberg form. 

Table 

3.8-5 

ATEICt 

Find eigenvalues of a matrix in 
Hessenberg form 

Table 

3.8-6 

ALPHA 

Compute integral of a function of 
the form (/ e^) * B * (/ eA) 

Table 

3.8-7 

MLINEQ 

Solve matrix equation A^x + A + C = 0 
for X 

Table 

3.8-8 

PRMR 

Print matrix routine 

Table 

3.8-9 

EQAT 

Set two matrices equal to each other 

Table 

3.8-10 

IDEN 

Set a matrix to the identity matrix 

Table 

3.8-11 

JUXC 

Concatenate two matrices columnwise, 
i . e . , A]b 

Table 

3.8-12 

JUXR 

Concatenate two matrices rowwise, 
i . e . , — 

Table 

3.8-13 

MNUS 

B 

Multiply matrix by (-1) 

Table 

3.8-14 

MULT 

Multiply two matrices 

Table 

3.8-15 

NORM 

Compute norm of a matrix 

Table 

3. 8-16 

SMLT 

Multiply a matrix by a scalar value 

Table 

3.8-17 

SUBT 

Subtract one matrix from another 

Table 

3.8-18 
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TABLE 3.8-1 

UTILITY ROUTINES (Continued) 


SUBROUTINE 

DESCRIPTION 

TABLE OF I/O 

TMLT 

Multiply a matrix by the transpose 

...... 


of another 

Table 3.8-19 

TRNP 

Compute transpose of a matrix 

Table 3.8-20 

ZERO 

Zeros a matrix 

Table 3.8-21 

SETA 

Compute ti'ansition matrix and dis- 
crete form of the noise matrix 

Table 3.8-22 

PROP 

Propagate covariance matrix 

Table 3.8-23 

EXPI 

Compute state transition and its- 
integral 

Table 3.8-24 

PRPA 

Propagate covariance matrix 

Table 3.8-25 

SJ.IPL 

Compute optimal Kalman gain matrix 

Table 3.8-26 

RCT 

and performs filter update 
Solve covariance equation 

Table 3.8-27 

INV 

Calculate inverse of a matrix 

Table 3.8-28 

RSET 

Compute transition matrix for con- 
tinuous Riccati equation 

Table 3.8-29 

EXPJ 

Compute exponential of a matrix 

Table 3.8-30 

ADD 

Add two matrices 

Table 3.8-31 

BETA 

Assist subroutine ALPHA in computing 
integral 

Table 3.8-32 

DINTEG 

Compute discrete integral 

Table 3.8-33 


TABLE 3.8-2 

INTERFACE TO ROUTINE DOT 


INPUTS 




N = 

dimension 

of 

vectors A and B 

A = 

N- vector 



B = 

N- vector 



OUTPUTS 



DOT = 

result of 

the 

dot product of A and B 
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TABLE 3.8-3 

INTERFACE TO ROUTINE CONTRL 

INPUTS 

N = number of states in system 
M = number of controls in system 
P = dimension of guidance command 
IP = run number 
C = MxN gain matrix 

OUTPUTS 

DELU = M- vector = (feedback gain)*(feedward control) 
DELV = P-vector of guidance inputs to the system 


TABLE 3.8-4 

INTERFACE TO ROUTINE STEADY 


INPUTS 

ITYP = flag to determine velocity or angular system 
PX = 201x31 time history computation matrix 
BELT = sample period 

IMAXl = number of time steps in sample periods, 
plus one 

OUTPUTS 

SS = steady-state value for variables 1 to NM3 
PV = peak value 
TPV = time of peak value 

T05 = time to reach a five percent bound of 
steady state 

T90 = time to reach ninety percent of steady state 


TABLE 3.8-5 

INTERFACE TO ROUTINE HSBG 

INPUTS 

N = order of the matrix 

A = the input matrix, N x N 
lA = size of the first dimension assigned to 
the array A in the calling program when 
the matrix is in double-subscripted data 
storage mode. IA=N wnen the matrix is in 
SSP vector storage mode. 

OUTPUTS 

A = returned N x N matrix in 'almost' Hessen- 
berg form 
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TABLE 3.8-6 

INTERFACE TO ROUTINE ATEIG 

INPUTS 

M = order of the matrix 
A = the input matrix, M by M 
lANA = vector whose dimension must be greater than 
or equal to M, containing on r«?turn indica- 
tions about the way the eigenvalues appeared 
( see Ref . 4 ) . 

lA - size of the first dimension assigned to the 

array A in the calling program when the matrix 
is in double-subscripted data storage mode. 
IA=M when the matrix is in SSP vector storage 
mode . 

OUTPUTS 

RR = vector containing the real parts of the eigen- 
values 

RI = vector containing the imaginary parts of the 
eigenvalues 


TABLE 3.8-7 

INTERFACE TO ROUTINE ALPHA 

INPUTS 

N = dimension of square matrix A 
A = constant matrix 
SIG = upper limit of integration 
EPSLO = convergence criterion 
GEST = index for starting test for convergence 
LLL =1,2 

H = 1.0, 0.0 

OUTPUTS 

SUM = the integral of the desired function 


TABLE 3.8-8 

INTERFACE TO ROUTINE MLINEQ 


INPUTS 

TOL = convergence tolerance 
N = dimension of square array 
A = NxN matrix 
C = NxN matrix 

OUTPUTS 

X = positive semi-definite solution of matrix 
linear equation 
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TABLE 3.8-9 

INTERFACE OF ROUTINE PRMR 


INPUTS 

A = Ilolerith string title of length L 
L = character length of Holerith field in 
N = mimber of rows 
M = number of columns 
X = NxM matrix 


TABLE 3.8-10 

INTERFACE TO ROUTINE EQAT 


INPUTS 


N = 

number of rows 

M = 

number of columns 

A = 

NxM matrix 

OUTPUTS 

B = 

NxM Matrix = A 


TABLE 3.8-11 

INTERFACE TO ROUTINE IDEN 


INPUTS 

N = dimension of matrix A 
OUTPUTS 

A = NxN identity matrix 


TABLE 3.8-12 

INTERFACE TO ROUTINE JUXC 


INPUTS 




N = 

number of rows in 

matrices A and B 

MA == 

number of colvunns 

in 

A 

MB = 

number of columns 

in 

B 

A = 

NxMA matrix 



B = 

NxMB matrix 



OUTPUTS 



C = 

Nx(MA + MB) matrix 
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TABLE 3.8-13 


INTERFACE TO ROUTINE JUXR 


INPUTS 
M = 

number of columns 

In matrices A and B 

NA = 

number of rows in 

A 

NB = 

number of rows in 

B 

A = 

NAxM matrix 


B = 

NBxM matrix 


OUTPUTS 


C = 

(NA+NB)xM ma,trix 



TABLE 3 . 8-14 

INTERFACE TO ROUTINE MNUS 
INPUTS 

N = number of rows 
M = number of columns 
A = NxM matrix 

OUTPUTS 

B = NxM matrix = -A 


TABLE 3.8-15 

^ INTERFACE TO ROUTINE MULT 


INPUTS 

N = number of rows in 
L = number of columns 
M = number of columns 
A = NxL matrix 
B = LxM matrix 

A 

A and number of rows in B 
in B 

OUTPUTS 

C = NxM matrix = A*B 



TABLE 3.8-16 

INTERFACE TO ROUTINE NORM 
INPUTS 

N = number of columns 
M number of rows 
A = NxM matrix 

OUTPUTS 

S = norm of matrix A 


fi,ffK0DUCD5ILrrY 0? THE 
OMCKNAL PAGE S' S ‘I'-r;, 








THE ANALYTIC SCIENCES CORPORATION 


TABLE 3,8-17 

INTERFACE TO ROUTINE SMLT 



INPUTS 

S = 

scalar to multiply A 

N = 

number of rows 

M = 

number of columns 

A = 

NxM matrix 

OUTPUT 

B = 

NxM matrix = S * A 


TABLE 3.8-18 

INTERFACE TO ROUTINE SUBT 



INPUTS 
N = 

number of rows 

M = 

number of columns 

A = 

NxM matrix 

B = 

NxM matrix 

OUTPUTS 

C = 

NxM matrix = A-B 


TABLE 3.8-19 

INTERFACE TO ROUTINE TMLT 



INPUTS 
N = 

number of rows in 

A 

L = 

number of columns 

in A and B 

M = 

number of rows in 

B 

A = 

NxL matrix 


B = 

MxL matrix 


OUTPUTS 


C = 

NxM matrix = A*B^ 



TABLE 3.8-20 

INTERFACE TO ROUTINE TRNP 



INPUTS 


N = 

number of rows 

M = 

number of columns 

A = 

NxM matrix 

OUTPUTS 

B = 

MxN matrix = A^ 
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TABLE 3.8-21 

INTERFACE TO ROUTINE ZERO 
INPUTS 

N = number of rows 

M = number of columns i 

OUTPUTS 

A = NxM zero matrix 


TABLE 3.8-22 

INTERFACE TO ROUTINE SETA 


INPUTS 

N - dimension of system 
F = NxN system matrix 
Q = NxN noise matrix 
DT = time step 

OUTPUTS 

RK = NxN discrete noise matrix 
TR = NxN transition matrix 


TABLE 3.8-23 

INTERFACE TO ROUTINE PROP 


INPUTS 

N = number of states 
P = NxN covariance matrix 
TR = (2*N) X (2*N) 

OUTPUTS 

P = NxN propagated covariance matrix 


TABLE 3.8-24 

INTERFACE TO ROUTINE EXP I 


INPUTS 

N = number of states 
F = NxN system matrix 
DT = specified time interval 

OUTPUTS 

TR = NxN transition matrix 
TRI = NxN integral matrix 
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TABLE 3.8-25 

INTERFACE TO ROUTINE PRPA 


INPUTS 

N = dimension of system 
P = NxN covariance matrix 
TR = NxN transition matrix 
QD = discrete noise matrix 

OUTPUT 

P = NxN propagated covariance matrix 


TABLE 3.8-26 

INTERFACE TO ROUTINE SMPL 


INPUTS 

N = dimension of covariance matrix 
M = dimension of measurement noise matrix 
P = NxN covariance matrix 
H = MxN measurement matrix 
R = measurement noise matrix 

OUTPUTS 

P = NxN updated covariance matrix 
GK = NxM optimal Kalman gain matrix 


TABLE 3.8-27 

INTERFACE TO ROUTINE RCT 


INPUTS 

N = dimension of covariance matrix P 
MAX = maximum number of time steps through 
which P is to be propagated 
P = NxN covariance matrix 
TR = (2+N) X (2*N) transition matrix 
EPS = convergence error criterion 

OUTPUTS 

P = NxN propagated covariance matrix 
NUM = number of time steps through which 
P has been propagated 
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TABLE 3.8-28 
INTERFACE TO ROUTINE INV 


INPUTS 

N = order of matrix A 

A = NxN matrix which is destroyed and replace 
by its inverse 

OUTPUTS 

A = NxN inverse of input matrix A 
D = Determinant of A -equals 0 if A is singular 




TABLE 3.8-29 

INTERFACE TO ROUTINE RSET 


INPUTS 

N = number of states 
M = dimension of measurement matrix 
F = NxN system matrix 
R = MxM m.easurement noise matrix 
H = MxN matrix relating measurement vector 
to state vector 
QG = NxN noise matrix 

OUTPUTS 

TR == (2*N) X (2*N) transition matrix 




TABLE 3.8-30 

INTERFACE TO ROUTINE EXPJ 


INPUTS 

N = dimension of matrix A 
A = NxN matrix 

OUTPUTS 

EXPA = NxN matrix exponential of A 


TABLE 3.8-31 
INTERFACE TO ROUTINE ADD 



INPUTS 
N = 

number of rows 

M = 

number of columns 

A = 

NxM matrix 

B = 

NxM matrix 

OUTPUTS 

C = 

NxM matrix = A+B 
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TABLE 3.8-32 

INTERFACE TO ROUTINE BETA 


INPUTS 


L = 

flag equal to 1 or 2 

N = 

size of matrix A 

A = 

NxN matrix 

ASIG = 

NxN matrix 

EPSLO = 

convergence criterion 

OUTPUTS 

SS = 

NxN matrix 


TABLE 3.8-33 

INTERFACE TO ROUTINE DINTEG 


INPUTS 

NT = number of summations 
N = dimension of system 
A = NxN matrix 
C = NxN matrix 

OUTPUTS 

S = NxN discrete integral matrix 
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4. SCHED USAGE 

SCHED transforms independent variables, performs 
polynomial or multiple regressions, and determines the 
correlation coefficients between gains, trim settings, and 
flight conditions. These features are controlled by the 
ICHPTR(I) flags (I = 1 to 6). The purpose of this program 
is to develop a set of regression coefficients to be used 
in an open-loop explicity-adaptive control algorithm. 


4 . 1 SCHED INPUT 

SCHED uses Fortran logical unit 5 (TAPES) for all 
its inputs. In creating the input deck for SCHED, it is 
necessary to choose appropriate flight conditions (inde- 
pendent variables X) and gains (dependent variables Y) which 
have been punched by DIGADAPT. Table 4.1-1 summarizes all 
card types for SCHED, including a list of variables, their 
descriptions, their format, and their appropriate card image 
positions. 
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TABLE 4.1-1 


INPUT FOR SCHED 


VARIABLE(S) 

DESCRIPTION 

FORMAT 

COLUMNS 

CARD 

TYPE 

N 

Number of observations 

14 

1-4 

1 

M 

Number of independent variables 

14 

5-8 

1 

L 

Number of dependent variables 

14 

9-12 

1 

ICHPTR 

ICHPTR(l) 

6-vector control flag 

=0, do not perform polynomial 

regression 

“1, performs polynomial regres- 
sion of highest order 2 of in- 
dependents and dependents 

12 

1-2 

2 

ICHPTR(2) 

=0, printout only correlation 
coefficients of polynomial 
regression 

•=1, print regression coeffi- 
cient, an analysis of variance, 
table of residuals and corre- 
lation coefficients 

12 

5-6 

2 

ICHPTR(3) 

=0, do not perform simple re- 
gression 

=1, perform simple regression 
between dependent variables 

12 

9-10 

2 

ICHPTR{4) 

=0, printout only correlation 
coefficients of simple regres- 
sion 

=1, print regression coeffi- 
cients, an analysis of variance, 
table of residuals, correlation 
coefficients, and an analysis 
of a multiple regression. 

12 

13-14 

2 

ICHPTR(5) 

=0, do not perform multiple 
regression 

=1, perform multiple regression 
analysis 

12 

17-18 

2 

ICHPTR(6) 

=0, printout only correlation 
coefficients of multiple re- 
gression 

=1, print regression coeffi- 
cients, analysis of multiple 
regression, table of residuals, 
and correlation coefficients. 

12 

21-22 

2 

NIND 

Number of independent variables 
for polynomial regression .< M 

12 

1-2 

3 

IIND 

NIND vector used to choose which 
independent variable will be 
used in polynomia.1 regressicin 

10(12, 2X) 

1-2, 5-6, 
9-10 , etc . 

4 

JORDER 

(NIND + 2) array of flags where 
the second subscript turns 'on' 
(=1) or 'off' (=0) the first 
and second-order polynomial re- 
gression. The first sub- 

script corresponds to the NIND 
variables of the IIND vector. 
These values are read in row- 
wise; hence, each row corre- 
sponds to an independent vari- 
able. 

2012 

1-2, 3-4, 
5-6 , etc. 

5 
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TABLE 4.1-1 
INPUT FOR SCHED 


VARIABLE(S) 

DESCRIPTION 

FORUAT 

COLUUNS 

CARD 

TYPE 

JIND 

(NIND X 5) array of flags which 
control the transformation on the 
independent variable. Each row 
of five switches corresponds to 
the appropriate independent 
variable of the IIND vector. 0 
Implies 'off and 1 Implies 'on'. 
The five switches pertain to the 
following transformations of tbs 
Independent variable 

3SI2 

1-2, 3-4 
5-6 , etc . 

6 


1. No change 

2 . Squared 

3. Inverse 

4. Inverse squared 

5. X • ABS(X) 




ISI 

Total number of multiple re- 
gressions to be performed 

12 

1-2 

7 

UUL 

Number of independent variables 
for each multiple regression. 
Card Type S and card type 9 is 
repeated ISI times in the se- 
quence 8, 9, 8, 9, etc. 

12 

1-2 

8 

UCC 

UUL - vector which specifies 
which independent variables 
are to be used for the multiple 
regression 

10(12,3X) 

1-2, 5-6 
9-10 , etc. 

9 

I HUB 

Case numbers (same as N'CASE 
from DIGADAPT). Cards type 

10, 11, 12 are repeated N times 
in the sequence 10, 11, 42, 10, 

11, 12, ... 

15 

1-5 

10 

X 

U independent variables 

6E12.S 

1-12,13-24 

11 

Y 

L dependent variables 

6E12.5 

1-12,12-24, 

12 


4.2 SCHED OUTPUT 

The output of SCHED is on Fortran logical unit 6 
(TAPE6). Table 4.2-1 gives a summary of the reports and a 
description and/or variable listing. All output from SCHED 
is printed from the main routine. 
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TABLE 4.2-1 


OUTPUT FROM SCHED 


REPORT 

DESCRIPTION-VARIABLE 

ProbJem/Parameter Cards 

N, number of observations, M, number of Independent 
variables and L, number of dependent variables 

Independent Variables for 
Polynomial Regressions 

IIND, JORDER, JIND 

Number of Multiple 
Regressions 

ISI 

Number of Independent 
Variables 

MUL, MCC 

List of Independent 

Case number and independent variables are printed 

Variables 

out along with headings for x, y, z, 6, i)i, v. 
velocity 

List of Dependent 
Variables 

Dependent variables y and the following variables: 

Cain Number 

I 

Mean 

AHEAN 

Standard Deviation 

STDV 

Percent 

PRCNT («100*STDV/AIIEAN) 

Begin Polynomial Re- 
gression 


Independent •= 

Independent variables modified 

Change = 

Change of variable code - IHNGE 

Ga In No 

Gain number - IAIN 

IRDER » 

Order of polynomial - IRDER 

Variable power 

I 

Moan 

XDAR 

Standard Deviation 

STD 

Correlation 

D 

Intercept 

ANS(l) 

Regression Coefficients 

B 

Multiple Correlation 
Cocf f iclont 

ANS(2) 

Source of Variation 

IRDER 

Degree of Freedom 

ANS(4) 

Sum of Squares 

ANS(6) 

Mean Square 

ANS(IO) 

Deviation About Re- 

NI , ANS(7) - array of square of deviation from 

gression 

regression, ANS(9) mean square of ANS(7) 

Total 

NT “ N-1, SUMSQ(NRDER) 

Table of Fiesiduals 

I, observation number 
XX(1), X value 
XX( IRDER » N + I), y value 
COE(NRDER), y estimate 
RES, residuals 

Independent Variable 

IIND(I) 

Polynomial Order 

NZ 

Correlation Coefficients 

CORRCO 

Begin Simple Regression 

Written on top of page 

Indejjondent and Dependent 
Variable 

XX 

Begin Multiple Regression 

Header Title 

Inttn-cept 

ANS(l) 

Multiple Correlation 

ANS(2) 

Standard Error of Estimate 

ANS(3) 

Multiple Regression No, 

INI 

Independent Varlubli = 

MCC, also prints out correlation coefficients CORMUL 
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5. 


DETAILED DESCRIPTION OF SCHED 


SCHED performs a correlation and regression 
analysis between independent and dependent variables. The 
independent variable.s are flight conditions at particular 
points in the flight envelope of the aircraft, and the de- 
pendent variables are control gains, trim settings, system 
matrices, and filter gains. 

In addition to performing regressions and deter- 
mining correlation coefficients, the program can perform 
transformations on any number of the independent variables. 

The change of variables that can be chosen are squared, 
inverted,, inverted and squared, and squared preserving the 
sign of the variable. The transformed variables are used 
to obtain a polynomial regression of order one or two. In 
addition, SCHED can perform simple regressions between the 
dependent variables and multiple regressions between a 
select group of independent variables and all the dependent 
variables. The outputs of the program are correlation 
coefficients, regression coefficients and tables of residuals. 
These results are computed by the following routines: GDATA, 

MULTR, CORRE, and ORDER. (See Fig. 5-1). 

GDATA generates independent variables up to the nth 
power (for the polynomial regression) and calculates means, 
standard deviations, sums of cross-products of deviations 
from means, and product-moment correlation coefficients (Ref.4). 
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n>i»7«i 






GDATA 


0»U generation 
for simple 
reg> ession 



Print regression 
coefficients and 
analysis of variance 



G> 


R-19766 


o 


CORBie 


Compute .means 
deviation and 
correlation 



Print table of 
residuals and 
correlation 
coefficients 


^ End MAIN ^ 


Figure 5-1 Flowchart of Program SCHED 
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For the polynomial regressions, ORDER chooses a de- 
pendent variable and a subset of independent variables, and 
MULTR computes the regression coefficients. For the multiple 
regressions, the routine CORRE is used in place of GDATA to 
calculate standard deviations, means, sums of cross-products 
of deviations from means, and product moment correlation coef- 
ficients from input data. 
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6 . 


TVHIS USAGE 


TVHIS perforins a flight path simulation based on a 
linear-tirae-varying model that includes actuator and rotor 
dynamics. The logic is controlled by the flags IREF, lAUG, 
lEST, lACT, IGNPT, and INOMPT. The purpose of this program 
is to test and evaluate the gain-scheduled control law table 
that has resulted from the analysis performed with DIGADAPT 
and SCHED. 

6.1 TVHIS INPUT 

Input for TVHIS is in two forms: card input on 

Fortran logical unit 5 (TAPES) and permanent file disk inputs 
on Fortran logical unit 1 (TAPEl). Table 6.1-1 describes 
the input variables, their descriptions, their units, their 
formats, and their appropriate card image positions. Refer- 
ence is made to Table 6.1-3 for the units of the units of 
the perturbation states. The variables ACONST and REGCOF 
should be obtained from the output of program SCHED. The 
disk file inputs, which are the same aerodynamic coefficients 
used in program DIGADAPT, are linearly interpolated in TVHIS 
to determine the system matrices FO and FI. 

In addition, there are variables initialized in 
BLOCK DATA. Table 6.1-2 gives the BLOCK DATA variables , 
their descriptions and their units. (Note that recompila- 
tion is necessary to change a variable in BLOCK DATA.) 
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TABLE 6.1-1 
CARD INPUT FOR TVHIS 


VARIABLES 

DESCRIPTION 

FORMAT 

UNITS 

COLUMN 

CARD 

TYPE 

XIC 

24-vector of Initial conditions 
on perturbation state vector 

16F4.0 

* 

1-4, 5-8, etc. 

1 

DIGUID 

Guidance time Interval 

F7.0 

sec 

1-7 

2 

DTCNTR 

Control time interval 

F7.0 

sec 

8-14 

2 

DTOUTP 

Output time interval 

F7.0 

sec 

15-21 

2 

NEV 

Number of special events 

13 

- 

1-3 

3 

ITYP.TEVENT 

EVUAG 

Key, time, and magnitude of 
special event - Card Type 4 is 
repeated NEV times 

I3,2F10.0 

-,sec,- 

1-3,4-7,8-11 

4 

NWPF.NWPL 

Flight waypoint numbers 

213 

- 

1-3, 4-6 

5 

ACONST 

68-vector of gain regression 
(9 cards) coefficients, from 
SCKED 

8F10.0 


1-10,11-20, 
21-30, etc. 

6 

REGCX3F 

68 X 3 array of u-regression 
coefficients, w-regression 
coefficients, and V^-regression 
coefficients from SCKED <26 
cards) 

8F10.0 


1-10,11-20, 
21-30, etc. 

7 


*See Table 6.1-3 


TABLE 6.1-2 

BLOCK DATA INITIALIZATION 


VARIABLES 

DESCRIPTION 

UNITS 

KIN 

Input logical unit = 5 

- 

KOUT 

Output logical unit = 6 

- 

KPUNCH 

Punch logical unit = 7 

- 

ACCU 

Convergence criterion for RKGS 

- 

NDIM 

Number of states = 24 

- 

WT 

NDIM vector of state error 
weights for RKGS 

- 

GVD 

Guidance gains at each waypoint 

sec”^ 

LIMIT 

Actuator rate and displacement 
limit key 

- 

GRAY 

Gravitational constant g 

4:4. / 2 

ft / sec 

RADS 

Conversion factor from degrees 
to radians 

rad/degree 

CKNOTS 

Conversion factor from knots 
to feet/sec 

feet/sec-knots 
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TABLE 6.1-2 

BLOCK DATA INITIALIZATION (Continued) 


VARIABLES 

■ — 

DESCRIPTION 

UNITS 

I XX 

Second moment of inertia about 
x-axis 

slug-ft^ 

lYY 

Second moment of inertia about 
y-axis 

slug-ft^ 

IZZ 

Second moment of inertia about 
z-axis 

slug-ft^ 

IXZ 

Product of inertia 

slug-ft^ 

XVECN 

Nominal state vector 

Same as state 
vector 

XIDTN 

Nominal inertial velocity 
vector 

ft/sec 

I REF 

=0 nominal state and control 
set to estimated nominal 
=1, nominal state and control 
set to true nominal 


I AUG 

=1, nominal state estimates 
used to get transformation 
matrix T 

=2, total state estimates used 
to get transformation 


lEST 

=0, for true nominal state 
used as independent regression 
variable 

=1, for the total state used 
as independent regression 
variable 


I ACT 

=0 actuator model off 
=1 actuator model on 

- 

UNCON 

4 vector static trim control 
regression intercepts 

stick equiva- 
lent inches 

UNREGC 

4x3 array static trim control 
regression coefficients 

- 

ANGCON 

2 vector static trim Euler 
angle regression intercepts 

degrees 

ANGRGC 

2x3 array static trim Euler 
angle regression coefficients 

- 
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TABLE 6.1-2 

BLOCK DATA INITIALIZATION (Continued) 


VARIABLES 

DESCRIPTION 

UNITS 

VIND 

Inertial velocity at each 
waypoint 

ft/sec 

KIND 

Inertial position at each 
waypoint 

feet 

TIMD 

Elapsed time at each waypoint 

sec 

IPLOT 

Not used 

- 

IGNPT 

=0, no gains printed out 
=1, gains printed out 

- 

INOMPT 

=0, do not print independent 
variable for regression 
=1, print independent variable 
for regression of nominal state 
and control values 



TABLE 6.1-3 
PERTURBATION STATES 


VARIABLES 

DESCRIPTION 

UNITS 

X 

Earth x-axis position 

feet 

y 

Earth y-axis position 

feet 

z 

Earth z-axis position 

feet 

u 

Body x-axis velocity 

ft/sec 

v 

Body y-axis velocity 

ft/sec 

w 

Body z-axis velocity 

ft /sec 

p 

Body x-axis angular rate 

rad/sec 

q 

Body y-axis angular rate 

rad/sec 

r 

Body z-axis angular rate 

ra d/sec. 

0 

Pitch Euler angle 

rad 


Roll Euler angle 

rad 
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TABLE 6.1-3 

PERTURBATION STATES (Continued) 


VARIABLES 

DESCRIPTION 

UNITS 

'>P 

Yaw Euler angle 

rad 


Differential collective (rotor) 

inches 

'x 

Collective (rotor) 

inches 

«s 

Lateral cyclic (rotor) 

inches 

^R 

Differential cyclic (rotor) 

inches 


Diff. collective rate (rotor) 

in/sec 

'6 

Collective rate (rotor) 

in/sec 

^S 

Cyclic rate (rotor) 

in/sec 


Differential cyclic rate (rotor) 

in/sec 

<Sba 

Diff, collective actuator 

inches 

'^Ca 

Collective actuacor 

inches 


Cyclic actuator 

inches 


Differential cyclic actuator 

inches 


6.2 TVH IS OUTPUT 

The output of TVHIS is stored on Fortran logical unit 6 
(TAPE6). Table 6.2-1 contains a summary of the output reports, 
a description and/or variable listing, and the subroutine 
which generates the report. 

TABLE 6.2-1 
OUTPUT FROM TVHIS 


REPORT 

DESCRIPTION- VARIABLE 

SUBROUTINE 

Integration Report 

TIMD - initial time 
TIMl - final time 
PRMT(3) - integration 
interval 

XVEC - state initial 
condition 

MAIN 
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TABLE 6.2-1 

OUTPUT FROM TVHIS (Continued) 


REPORT 

DESCRI PT I ON- VARI ABLE 

SUBROUTINE 

Input Report 

DTGUID - guidance time in- 
terval 

DTCNTR - control time in- 
terval 

DTOUTP - output time inter- 
val 

RDDAT 

Special Event 
Specification 

ITYP - key to special event 
TEVENT - time of special 
event 

EVMAG - magnitude of spe- 
cial event 

RDDAT 

Flight Waypoint No. 

NWPF, NWPL 

RDDAT 

Gain Regression 
Coefficients 

ACONST, REGCOF 

RDDAT 

Initial Waypoint 
System Matrix 

FO 

MISION 

Initial Waypoint 
Input Matrix 

GO 

MISION 

Initial Waypoint 
Nominal Inertial 
Velocity 

XIDTO 

MISION 

Initial Waypoint 
Nominal State 

XVECO 

MISION 

Final Waypoint 
System Matrix 

FI 

MISION 

Final Waypoint 
Input Matrix 

G1 

MISION 

Rows = 
Column = 

Print matrix routine out- 
put description head plus 
number of row and columns 

PRMR 

Special Event 
Report 

TIME 

OUTP 

DYDS = 

Guidance Off 

Guidance perturbation 
input 

OUTP 

DXS = 

Guidance Off 

State perturbation input 

OUTP 
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TABLE 6.2-1 

OUTPUT FROM TVHIS (Continued) 


REPORT 

DESCRIPTION- VARIABLE 

SUBROUTINE 

DVS = 

Guidance Off, 
Controller Off 

Control perturbation input 

OUTP 

Time Histories 

The total state and control 
vectors (PX) are x, y, z, u, 
V, w, p, q, r , <t>, .ip, 6b, 

6C, 6s, 6r, 6b, 6c, js > ^R> 

*^B>j' > *^R^> ^B<p > 6Crji, 

u, V, w,xi,yj, Zj. 

OUTP 


Similarly, perturbation state 
and control vectors (PDX) 
are output 


Control Report 

TIME - time in seconds 
V - augmented error state 
UCOM - control command 

CNTRL 

State Feedback 
Gains 

CKX 

CNTRL 

Integrated Feedback 
Gains 

CKV 

CNTRL 

Control Feedback 
Gains 

CKU 

CNTRL 

X NOM EST 
U NOM EST 

XVNEST - X nominal esti- 
mate 

CNTRL 

Perturbation States 
and Control 

DX, DU 

CNTRL 

Control For Appli- 
cation 

UCOM 


Integrated Output 
States 

GDCMVC - guidance comniahd 
vector 

V - augmented error state 
DUM2 - approximate time 
derivative of the aug- 
mented error state 

CNTRL 

Transformation 
Matrix T 

T - transformation from state 
vector to output vector 

TCO-MP 

Guidance Report 

TIME 

GUIDE 

GDCMVC 

Guidance command vector 

GUIDE 

Turbulence Included 

Variation from subroutine 
TURC included 

TURC 
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7 . DETAILED DESCRIPTION OF TVHIS 


TVHIS is designed as a realistic yet efficient 
helicopter simulation, it is based on a linear-time- 
varying helicopter model, including actuator and rotor 
dynamics models. The capability for simulating rate and 
displacement limits on the actuators also is incorporated. 
TVHIS includes a guidance subroutine which calculates total 
inertial velocity commands by summing the nominal velocity 
(specified by a nominal trajectory) and a velocity correc- 
tion due to cross range and altitude errors. The guidance 
algorithm may be bypassed to allow testing of the control 
algorithm alone. TVHIS includes a digital-adaptive control 
routine that includes a state estimation section, a gain 
adaptation procedure, and the controller logic. The compu- 
tations in TVHIS are accomplished by the following major 
subroutines: RDDAT, MIS ION, FLTINT, FGCOMP, RKGS , FCT , OUTP , 

CNTRL and GUIDE (see Figure 7-1). 

RDDAT is the input routine which reads the flags 
and input parameters and which prints them out. The purpose 
of MISION is to pick the waypoints of the trajectory where 
position and velocity have been specified. Using the flight 
conditions, FLTINT linearly interpolates between the aero- 
dynamic coefficients and then FGCOMP computes the system 
matrix. 

RKGS has been incorporated from the IBM Scientific 
Subroutine Package (Ref. 4). The purpose of the Runge-Kutta 
method is to obtain an approximate solution of a system of 
first-order ordinary differential equations with given 
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Figure 7-1 Overview of Program TVHIS 


initial values. It is a fourth-order integration procedure 
which is stable and self-starting; that is, only the func- 
tional values at a single previous point are required to 
obtain the functional values ahead. For this reason, it is 
easy to change the step size (h) at any step in the calcula- 
tions. Control of accuracy and adjustment of h is done 
by comparison of the results due to double and single step- 
sizes, 2h and h, 

FCT evaluates the differential equation using the 
state and control vectors. OUTP is used to apply special 
events, to produce time history tables, and to call guidance 
and control algorithms, which are found in subroutine GUIDE 
and CNTRL, respectively. 


Ei5?M3DUClBILI73( OF 
OBiGimAL - 
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GUIDE computes the vehicle guidance commands by 
combining the nominal command with a perturbation command 
due to inertial position error. The guidance algorithm 
may use time-varying gains. 

CNTRL includes the state estimation logic and the 
controller algorithm. This logic computes the perturbation 
states, which are the differences between the nominal and 
total state estimates. If desired, the true stste values 
may be taken as the estimates. The controller algorithm 
updates the control gains by applying the results of the con- 
troller gain regression of SCHED and calculates the control 
commands by applying the feedback law. 
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APPENDIX A 

IMPLEMENTATION OF PERMANENT FILE FOR 
AERODYNAMIC COEFFICIENTS 


The program which stores aerodynamic coefficients 
(AERO) must be run before DIGADAPT and TVHIS can be used. 

The purpose of this program is to create a permanent file 
on a disk. The file consists of 297 different cases of aero- 
dynamic coefficients whose parameter variations are explained 
below . 

The aerodynamic coefficients, corresponding to 

• • • 

specific values of ip, x and z stored on the permanent disk 
file, are selected by computing a case number (NCASE) 
defined as 

NCASE = (J-D+99 + (K-l)*ll + L 

where J is an index of 3 values for ip , K is an index of 9 
values for z, and L is an index of 11 values for x. The ^ 
(radian/sec) values of 0.0, +0.05, -0.05 correspond to 
J =1 to 3. The z (feet /min) value of -2000, -1500, -1000, 
-500, 0, 500, 1000, 1500, 2000 correspond to K=1 to 9. The 
X (knots) values of -40, -20, 0, 20, 40, 60, 80, 100, 120, 

140, 160 correspond to L = 1 to 11. For ea.ch case, there are 
71 aerodynamic coefficients used in the formulation of the 
system matrix F. 
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APPENDIX B 

FLOWCHARTING SYMBOLS 

Table B-1 gives the flowcharting symbols used in 
this manual. They have been defined in the ANSI standard 
(Ref. 5 and 6), with the exception of the off-page connec- 
tor, which is an IBM defined symbol, and the loop symbol, 
which is a TASC defined symbol. 

The symbol used within the off page connector, as 
in Fig. 3.1-1, identify the connector (by letter) and the 
connecting segment (by number in parentheses). 

TABLE B-1 
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APPENDIX C 
DESIGN PROCEDURES 


The design procedure for determining a control law 
using the programs DIGADAPT and SCHED is outlined below 
(Ref. i). The procedure indicates at each step which method 
was chosen for this report; alternate steps which could have 
been attempted are also provided. 


• Sampling Interval 

Decide on covariance matrix for gust states. 
Add any other system noises. 

Specify the state vector bounds. 

Propagate state covariance for various flight 
conditions until the state vector bound is 
reached. 

• Flight Conditions for Point Design 

Choose a representative number of the most 
common flight conditions. 

Separate the flight conditions into two 
groups, one near hover, the other at higher 
speeds. Near hover, constant gains can be 
chosen; mode switching occurs to the regressed 
gains after a specified forward velocity 
is exceeded. 

• Command-Controller 

Choose the command system and controller 
type for continuous or discrete case, e.g.. 
Velocity Command (PII Controller), 

Attitude Command (PI Controller), or 
Dynamic Controller (Regulator). 

• Criteria 

Specify suitable step response requirements, 
closed-loop eigenvalues, or 
quadratic weightings. 
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Determine the state variables to be 
weighted by quadratic synthesis. 

Specify Q and R adjustment policy if 
required. 

• Controller Gain Scheduling 

Obtain control gains from DIGADAPT. 

Iterate Q,R as required. 

Punch out controller gains and flight 
conditions . 

Choosf flight conditions for regression 
(body axes or inertial axes). Reorder 
punched cards with chosen flight condition 
type . 

Run SCHED for all combinations of independent 
variables (flight conditions), printing out 
correlation coefficients. 

Using SCHED output choose independent 
variables for regression based on correlation 
coefficients and availability of measurements 
onboard the vehicle. 

Run SCHED for regressions with chosen inde- 
pendent variables. 

Record regression coefficients for use in 
open-loop explicit-adaptive gain scheduling 
algorithm. 

• Trim Values 

Run DIGADAPT for as many flight conditions 
as desired, punching out control trim posi- 
tions and state positions. 

Choose body or inertial axes for Independent 
variable conditions. 

Schedule static trim values against flight 
condition . 

Schedule dynamic trim coefficients against 
flight condition. 

• Low-Pass Filters 

Establish signal and noise variances for 
measurement states. 

Define filter gains using discrete-time 
Kalman filter algorithms. 
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• Complementary Filters . 

Choose a sampling time for the filters if 
they are to be discrete. Filters are simple 
enough to be adjusted on line or by full 
simulation techniques. 

m Kalman Filter Design for Unmeasured States 

Define states which are measured, disturbance 
covariance matrix, and observation noise 
covariance matrix. 

Run DIGADAPT to obtain Kalman filter gains 
for flight conditions desired. Run SCHED with 
Kalman gains as noted above. 
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APPENDIX D 

LIST OF DELIVERABLES 


Delivered with this document are thirteen distinct, 
computer-related items: one magnetic tape, six executable 

card decks, and six computer-printed listings. These thirteen 
items are detailed in the following list. 


1. Magnetic tape, 7-track, 556 bpi (HI) containing 
five files. 

File 1 - (binary) DIGADAPT program library 
File 2 - (binary) SCHED program library 
File 3 - (binary) TVHIS program library 
File 4 - (binary) AERO program library 
File 5 - (coded) aerodynamic coefficient data 


2. Executable deck which reads item 1 (above) and 
creates four permanent files. 

DIGOBJ - DIGADAPT object code 
SCHOBJ - SCHED object code 
TVHOBJ - TVHIS object code 

AERO - data file of aerodynamic coefficients, 
read by DIGADAPT and TVHIS 


3. 

Executable 

deck 

- DIGADAPT sample 

execution 

#1 

4. 

Executable 

deck 

- DIGADAPT sample 

execution 

#2 

5. 

Executable 

deck 

- DIGADAPT sample 

execution 

#3 

6. 

Executable 

deck 

- SCHED sample execution 


7. 

Executable 

deck 

- TVHIS sample execution 
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8. Computer printout generated by execution of item 
2 - includes compiler-generated source listings 
of DIGADAPT, SCHED, TVHIS, and AERO and the 
execution- time output of program AERO. 


9. Computer printout generated by execution of item 3 - 
DIGADAPT sample output #1. 


10. Computer printout generated by execution of item 4 - 
DIGADAPT sample output #2. 


11. Computer printout generated by execution of item 5 - 
DIGADAPT sample output #3. 


12. Computer printout generated by execution of item 6 - 
SCHED sample output . 


13. Computer printout generated by execution of item 7 - 
TVHIS sample output. 


D-2 


THE ANALYTIC SCIENCES CaRPORATION 


REFERENCES 


1. Stengel, R.F, , Broussard, J.R. and Berry, P.W., 

"The Design of Digital Adaptive Controllers for VTOL 
Aircraft," The Analytic Sciences Corporation, 

TR-640-1, Reading, Massachusetts, September 1975. 

2. CYBERNET Scope 3.3 Reference Manual - Batch and Remote 
Batch Operating System, Publication No. 8661800, CDC 
CYBERNET Publications Department, Minneapolis, 

Minnesota, 1973. 

3. Sandell, N.R. and Athans, M. , "Modern Control Theory: 
Computer Manual," M.I.T, Center for Advanced Study, 
Cambridge, Massachusetts, 1974. 

4. System/360 Scientific Subroutine Package, Version III, 
Publication No. H20-0166-5, IBM Technical Publications 
Department, White Plains, New York, March 1970, 

5. American National Standard Flowchart Symbols and Their 
Usage in Information Processing, Publication No. 
AIM1171/375, ANSI, Inc., New York, New York, 1971. 

6. IBM Flowcharting Techniques, Publication No. C20-8152-1, 
IBM Technical Publications Department, White Plains, 

New York, March 1970. 



